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ABSTRACT 


Wireless  communication  technology  has  become  a  critical  aspect  in  many  civilian  and 
military  applications.  With  regard  to  remote  sensing,  search  and  rescue,  disaster  relief 
operations  and  signals  intelligence,  there  exists  an  interest  in  developing  capabilities  to 
collect  these  signals-of-interest.  The  objective  of  this  dissertation  is  to  maximize  signal 
collection  performance  in  the  presence  of  signal  measurement  and  sensor  related  errors. 
To  accomplish  this  objective,  we  proposed  a  signal  collection  scheme  that  exploits  an 
elevated,  mobile  network  to  maximize  the  collaborative  collection  of  a  target  signal. 

The  proposed  scheme  begins  with  source  localization.  This  technique  consists  of 
an  initial  weighted  least-squares  estimate  followed  by  a  maximum-likelihood  estimate. 
Implemented  on  an  elevated,  mobile  network,  this  technique  is  able  to  obtain  an  optimal 
localization.  To  enhance  localization  robustness,  we  developed  an  outlier  rejection 
process  that  mitigates  the  effects  of  measurement  and  sensor  position  errors. 

To  collect  the  signal,  this  research  quantified  the  effects  of  sensor  position  errors 
on  beamforming  and  proposed  a  novel  signal  collection  scheme  that  combines  signal 
estimation  and  collaborative  beamfonning.  Using  all  these  techniques  in  concert,  we 
were  able  to  show  that  the  proposed  scheme  outperforms  standard  collaborative 
beamforming  in  the  presence  of  sensor  position  errors. 
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EXECUTIVE  SUMMARY 


Wireless  communication  technology  has  penetrated  many  aspects  of  both  civilian  and 
military  applications.  With  cellular  technology  constantly  evolving  and  proliferating,  we 
now  see  wireless  communications  as  an  integral  part  of  daily  life.  With  regard  to  remote 
sensing,  search  and  rescue,  disaster  relief  operations  and  signals  intelligence,  there  exists 
an  interest  in  developing  capabilities  to  collect  these  signals.  The  objective  of  this 
dissertation  is  to  maximize  signal  collection  perfonnance  in  the  presence  of  various 
signal-  and  sensor-related  errors.  To  accomplish  this  objective,  we  proposed  a  signal 
collection  scheme  that  exploits  an  elevated,  mobile  network  to  maximize  the 
collaborative  collection  of  a  target  signal. 

With  its  roots  dating  back  to  the  development  of  radar  during  World  War  II, 
beamforming  has  found  a  prominent  place  in  modern  wireless  communications.  In  its 
classical  use,  through  the  manipulation  of  array  weights,  beamforming  coherently 
amplifies  a  signal  in  a  given  direction,  while  reducing  undesirable  signals  in  others. 
Although,  beamforming  is  a  powerful  technique  that  enables  spatial  filtering,  it  has  a 
strong  dependence  on  the  a  priori  knowledge  of  the  target  signal’s  direction  of  arrival 
[1].  Because  of  this  dependency,  beamforming  is  often  initialized  by  source  localization. 

Our  proposed  scheme  begins  with  a  source  localization  technique;  this  technique 
is  used  to  determine  the  signal’s  position  and  direction  of  arrival.  This  knowledge  is  then 
used  by  a  collaborative  beamformer  to  maximize  and  collect  the  target  signal.  Both  of 
these  techniques  are  enhanced  through  the  use  of  an  elevated,  mobile  network.  This 
mobility  allows  for  the  reconfiguration  of  the  sensor  network’s  topology  to  create  an 
ideal  sensor-target  geometry.  This  geometry  minimizes  any  geometric  dilutions  of 
precision,  thus  allowing  for  an  optimal  location  estimate  to  be  achieved  [2].  The 
network’s  elevation  increases  signal  power  and  range  while  enabling  unique  sensor 
formations  for  signal  collection  that  is  robust  against  sensor  position  errors. 

The  proposed  scheme  begins  with  a  localization  technique.  In  this  localization, 
we  propose  a  two-stage  technique  capable  of  obtaining  an  accurate  localization  in  the 
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presence  of  signal  noise  and  sensor  position  errors.  This  technique  consists  of  an  initial 
weighted  least-squares  estimate  followed  by  a  refining  maximum-likelihood  estimate. 
Here,  the  initial  estimate  is  used  to  reconfigure  the  network’s  topology.  This 
reconfiguration  is  used  to  create  an  optimal  sensor  network  to  target  geometry,  thus 
creating  the  necessary  conditions  for  the  refining  estimation  to  deliver  an  optimal  location 
estimate.  Our  simulation  results  show  the  proposed  localization  technique  to  be  efficient, 
i.e.,  it  approaches  the  Cramer-Rao  lower  bound  in  the  small  error  region  [2], 

To  enhance  the  localization  performance,  we  developed  a  measurement  outlier 
rejection  process  to  mitigate  the  effects  of  measurement  and  sensor  position  errors.  This 
technique  uses  a  combination  of  single  case  diagnostics  and  the  Mahalanobis  distance  to 
identify  and  remove  specious  measurements  for  the  least-squares  estimate  [3].  Through 
simulation,  we  demonstrated  this  technique  to  be  effective  in  the  presence  of  both 
measurement  and  position  errors. 

To  develop  the  proposed  signal  collection  scheme,  we  began  by  analyzing  the 
effects  of  sensor  position  errors  on  the  array  factor.  We  derived  an  expression  for  the 
effects  of  position  error  on  the  array’s  main  beam  gain.  Through  simulation,  we 
validated  our  results  and  showed  that  the  mean  value  of  the  main  beam  signal  phase  was 
unaffected  by  position  errors.  To  enhance  the  collection  performance  in  the  presence  of 
these  errors,  we  developed  a  signal  estimator  for  both  uniform  and  Gaussian  position 
errors.  Our  simulation  results  showed  a  maximum  improvement  in  array  gain  of 
approximately  37  percent  for  the  standard  deviation  of  position  error  values  greater  than 
0.4  m.  Other  results  demonstrated  that  our  approach  of  combining  a  signal  estimation 
technique  with  collaborative  beamforming  is  a  viable  means  of  collecting  a  target  signal 
from  an  elevated,  mobile  wireless  sensor  network. 

Taking  advantage  of  the  concept  of  a  unique  elevated,  mobile  wireless  sensor 

network  realized  through  the  use  of  multirotor  UAVs,  our  scheme  used  two  existing 

localization  techniques  to  deliver  a  precise  source  location  estimate.  Using  this 

information  with  statistical  knowledge  of  the  sensor  position  errors  in  a  signal  estimator, 

we  were  able  to  derive  a  novel  collaborative  signal  collection  scheme.  This  scheme  was 

shown  to  be  capable  of  collecting  and  amplifying  a  target  signal  in  the  presence  of  such 
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errors.  With  all  these  techniques  in  concert,  the  objective  of  this  dissertation  to  maximize 
signal  collection  performance  in  the  presence  of  various  signal-  and  sensor-related  errors 
was  achieved. 
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I.  INTRODUCTION 


Wireless  communication  technology  has  penetrated  many  aspects  of  both  civilian 
and  military  applications.  The  International  Telecommunications  Union  (ITU)  estimates 
that  there  are  approximately  seven  billion  cellular  subscribers  worldwide  [1].  This 
translates  to  a  global  penetration  rate  of  95.5  percent  and  is  expected  to  increase  each 
year.  Clearly,  wireless  communications  has  become  a  common  aspect  of  human  life. 
With  regard  to  remote  sensing,  search  and  rescue,  disaster  relief  operations  and/or  signals 
intelligence,  there  exists  an  interest  in  developing  capabilities  to  collect  these  signals. 
One  such  solution  is  collaborative  beamforming  [2], 

With  its  roots  dating  back  to  the  development  of  radar  during  World  War  II, 
beamforming  has  found  a  prominent  place  in  modern  wireless  communications. 
Beamforming  is  much  akin  to  the  matched  filter  but  from  a  spatial  perspective  [3]. 
Instead  of  linearly  weighting  a  time  domain  signal  for  signal  filtering,  a  beamfonner 
linearly  weights  spatial  array  data  for  spatial  filtering.  In  its  classical  use,  beamforming 
coherently  amplifies  a  signal  in  a  given  direction  while  reducing  undesirable  signals  in 
other  directions.  Although  beamforming  is  a  powerful  technique  that  enables  spatial 
filtering,  it  has  a  strong  dependence  on  a  priori  knowledge  of  the  target  signal’s  direction 
of  arrival  [4],  Because  of  this  sensitivity,  beamforming  is  often  preceded  by  a  direction 
of  arrival  estimation  phase. 

With  regard  to  direction-of-arrival  estimates,  hyperbolic  localization  [5]  has 
appealing  qualities  of  interest  to  this  research.  Of  particular  interest  is  hyperbolic 
localization  perfonned  from  an  elevated,  mobile  wireless  sensor  network  (WSN).  The 
tenn  hyperbolic  is  used  to  denote  the  use  of  signal  time-difference-of-arrival 
measurement  and  calculations.  The  network  is  elevated  to  promote  a  line-of-sight 
communication  link  for  better  range  and  reliability,  while  the  mobility  allows  for  the 
reconfiguration  of  sensor-target  geometries  for  more  precise  direction-of-arrival 
estimates. 
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The  extension  of  beamfonning  to  wireless  sensor  networks  has  given  rise  to  the 
concept  of  collaborative  beamforming.  Here,  a  fully  synchronized  network  of  sensors  is 
used  to  coherently  collect  and  combine  multiple  received  signals.  These  signals  along 
with  network  topology  infonnation  are  used  to  form  a  distributed  beamforming  array. 
This  technique  coupled  with  recent  advances  in  micro  sensor  and  multirotor  unmanned 
aerial  vehicle  (UAV)  technologies  has  led  to  the  concept  of  a  semi-stationary  airborne 
sensor  network.  Such  a  network  is  capable  of  providing  both  the  elevation  and  mobility 
needed  for  highly  accurate  hyperbolic  localization  while  providing  a  suitable  network 
platform  for  collaborative  signal  collection. 

A.  OBJECTIVE 

The  dissertation  objective  is  to  maximize  signal  collection  performance  in  the 
presence  of  various  signal  and  sensor  related  errors.  To  accomplish  this  task,  we  propose 
a  new  scheme  for  signal  collection  that  is  preceded  by  a  localization  phase. 

To  improve  the  scheme’s  localization  accuracy,  we  propose  the  use  of  an 
optimally  positioned  two-stage  localization  technique.  This  technique  consists  of  an 
initial  weighted  least-squares  estimate  followed  by  a  refined  maximum-likelihood 
estimate.  Here,  the  initial  estimate  is  used  to  reconfigure  the  network’s  topology,  thereby 
creating  an  optimally  configured  network  to  obtain  the  refined  hyperbolic  estimate.  To 
combat  the  effects  of  noise,  we  propose  the  use  of  measurement  outlier  diagnostics  in 
concert  with  the  two-stage  technique.  Here,  we  postulate  that  the  identification  and 
removal  of  specious  time-difference-of-arrival  measurements  will  improve  the  robustness 
of  localization  in  the  presence  of  noise. 

To  maximize  the  signal  collection  performance  in  the  presence  of  noise,  we 
propose  the  use  of  signal  estimation  in  conjunction  with  collaborative  beamforming. 
Using  the  concept  of  sensor  stacks  with  knowledge  of  the  sensor  positional  error 
statistics,  we  can  derive  a  signal  estimator  to  combat  the  effects  of  sensor  position  errors. 
To  further  enhance  the  collection  capability  for  noise  and  interference  rejection,  we  also 
propose  the  use  of  virtual  filling  to  increase  the  array’s  ability  to  reject  interfering  signals. 
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B. 


RELATED  WORK 


Beamforming  has  been  shown  to  be  an  effective  method  for  signal  collection  and 
interference  rejection  [6],  but  it  has  been  shown  to  be  highly  susceptible  to  array  steering 
vector  errors  [7],  This  is  when  an  array  steering  vector  is  not  in  line  with  that  of  the 
target  signal.  To  make  beamforming  more  robust  against  these  array  mismatch  errors, 
Ahmed  and  Evans  [8]  suggest  the  use  of  inequality  constraints  on  the  array  weights  while 
Er  and  Cantoni  [9]  suggest  the  use  of  derivative  constraints.  Both  of  these  approaches, 
although  effective,  also  decrease  the  array’s  interference  rejection  capabilities.  Li,  Stoica 
and  Wang  [10]  developed  a  robust  beamformer  that  minimized  the  effects  of  non-random 
steering  errors  using  a  diagonal  loading  method.  Lee  and  Lee  [11]  proposed  a  robust 
beamformer  for  signal  collection,  which  minimizes  a  cost  function  based  on  received 
signal  data  and  knowledge  of  steering  error  statistics.  The  signal  collection  scheme 
proposed  in  this  research  uses  elevated,  mobile  sensor  network  formations  and  signal 
estimation  in  the  context  of  collaborative  beamforming.  Preceded  by  a  localization  phase 
that  creates  a  unique  sensor  network-to-target  formation,  the  scheme  is  able  to  isolate  the 
effects  of  phase  perturbations  due  to  position  errors.  It  then  compensates  for  them 
through  the  use  of  sensor  formations,  signal  estimation,  and  knowledge  of  sensor  position 
error  statistics.  Unlike  the  schemes  of  Ahmed  and  Evans  [8]  and  Er  and  Cantoni  [9],  our 
scheme  does  not  require  any  constraints  on  the  beamformer  weights.  Similar  to  that  of 
Lee  and  Lee  [11],  our  scheme  utilizes  sampled  array  data  and  knowledge  of  error 
statistics  to  compensate  for  array  phase  errors.  The  difference  being  that  our  scheme  uses 
the  sampled  data  in  a  signal  estimation  context  with  consideration  to  array  formation 
rather  than  an  array  weight  optimization  problem  with  no  consideration  to  array 
formation. 

One  of  the  earliest  hyperbolic  location  systems  was  investigated  by  Carson  [12]  in 
1972.  He  proposed  a  system  of  hydrophones  to  collect  signal  arrival  times,  with  a  goal  to 
estimate  the  location  of  subsurface  acoustic  emitters.  Today,  numerous  extensions  and 
breakthroughs  have  been  made  to  adapt  these  concepts  for  radio  frequency  (RE)  signals. 
Pang  [13]  provides  a  closed  form  solution,  but  the  estimator  is  not  efficient,  i.e.,  it  does 
not  achieve  the  Cramer-Rao  lower  bound  [14],  [15],  Hahn  and  Tretter  [16]  used  the 
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hyperbolic  equations  in  a  maximum-likelihood  estimation  scheme  that  is  unbiased  and 
asymptotically  efficient.  Being  a  linearized  system  of  non-linear  equations,  the 
maximum-likelihood  estimator  requires  an  initial  reference  estimate  to  avoid  local 
minima.  Chan  and  Ho  [5]  propose  a  weighted  least-squares  approach  that  yields  an 
efficient  estimator  in  the  small  error  region.  Our  approach  will  combine  the  maximum- 
likelihood  estimate  of  Hahn  and  Tretter  [16]  and  the  weighted  least-squares  estimate  of 
Chan  and  Ho  [5]  to  create  a  two-stage  localization  technique.  Furthermore,  this 
technique  will  use  optimal  sensor  network  formations  as  outlined  by  Ho  and  Vicente 
[17].  This  combination  of  hyperbolic  estimation  and  optimal  sensor  network  formations 
will  be  used  to  provide  highly  accurate  localization  from  an  elevated,  mobile  WSN 
platfonn. 

Rousseeuw  and  Leroy  [18]  provided  the  first  comprehensive  review  of  robust 
estimation  and  the  effects  of  statistical  outlier  measurements.  Key  contributions  by 
Mahalanobis  [19]  and  Cook  [20]  helped  quantify  the  effects  of  outliers  in  estimation  by 
defining  the  Mahalanobis  and  Cook  distance.  Picard  and  Weiss  [21]  coupled  these 
findings  with  the  concept  of  sparse  representation  to  derive  a  bound  on  the  number  of 
outliers  that  can  be  identified  in  a  hyperbolic  localization  problem.  Many  researchers 
have  contributed  techniques  to  combat  the  effects  of  outlier  measurements.  McGuire, 
Plataniotis  and  Venetsanopoulos  [22]  proposed  the  use  of  data  fusion  techniques.  Yang, 
Wang  and  Luo  [23]  showed  that  convex  relaxation  methods  can  be  used  to  overcome 
outliers  induced  by  both  measurement  noise  and  sensor  position  errors.  Choi,  et  al.  [24] 
used  recursive  filtering  to  deal  with  errors  from  incorrect  stochastic  infonnation  in 
location  estimators.  With  the  proposed  localization  technique  implemented  on  an 
elevated,  mobile  sensor  network,  there  will  be  increased  measurement  noise  due  to  sensor 
position  errors.  To  minimize  the  effects  of  these  errors  as  well  as  measurement  outliers, 
we  propose  an  outlier  rejection  process  based  on  single  case  diagnostics  [18]  and  the 
squared-Mahalanobis  distance  [19].  Overall,  a  robust  localization  scheme  that  combines 
the  proposed  two-stage  localization  technique  and  the  outlier  rejection  process  is 
achieved. 
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C.  ORGANIZATION 

This  dissertation  is  organized  as  follows. 

The  fundamental  principles  involved  with  hyperbolic  localization  and 
beamforming  is  the  focus  of  Chapter  II.  Localization  based  on  the  time-difference-of- 
arrival  measurements,  weighted  least-squares  estimator,  maximum-likelihood  estimator, 
and  robust  signal  outlier  detection  techniques  are  discussed  in  the  localization  section. 
With  collaborative  signal  collection  as  the  end  goal,  we  also  provide  an  introduction  to 
beamforming  and  its  extension  to  collaborative  beamforming. 

The  proposed  two-phase  scheme  is  presented  in  Chapter  III.  This  includes  a 
description  of  the  scheme’s  process  flow  as  it  moves  through  the  two  phases  of 
localization  and  collection.  A  discussion  of  the  deployment  and  operating  concepts  is 
provided  to  give  context  to  the  proposed  scheme.  The  many  figures  of  merit  used  to 
evaluate  the  scheme’s  performance  are  also  detailed. 

The  detailed  analysis  in  support  of  the  robust  localization  phase  of  our  scheme  is 
included  in  Chapter  IV.  Additional  robust  signal  processing  techniques  that  are 
implemented  in  both  estimators  to  combat  the  effects  of  measurement  outliers  and  sensor 
position  errors  are  also  included.  The  performance  of  the  robust  localization  phase  is 
then  examined  using  simulation. 

The  signal  collection  phase  of  the  proposed  scheme  is  described  in  Chapter  V. 
This  includes  an  analysis  of  sensor  position  errors  and  their  effects  on  an  array’s  array 
factor.  Based  on  these  findings,  we  derive  a  novel  signal  estimator  to  combat  the  effects 
of  uniform  and  Gaussian  sensor  position  errors.  The  perfonnance  of  the  resulting 
collection  scheme  is  then  investigated  using  simulation. 

In  conclusion,  a  summary  of  the  dissertation,  a  discussion  of  the  research’s 
significant  contributions,  and  a  list  of  potential  future  research  topics  are  given  in  Chapter 
VI.  An  appendix  is  also  included,  which  contains  selected  MATLAB  simulation  scripts. 


5 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


6 


II.  BACKGROUND 


To  perform  signal  collection,  we  must  understand  the  different  processes  involved 
in  source  localization  and  collection.  In  this  section,  we  provide  the  necessary 
background  material  to  properly  present  our  research  and  results.  We  begin  with  a  brief 
discussion  on  WSN  that  forms  the  basis  of  the  proposed  signal  collection  scheme.  With 
the  premise  of  sensor  networks  explored,  we  discuss  the  principles  of  hyperbolic/TDOA 
source  localization  and  collaborative  beamforming. 

Throughout  the  dissertation,  the  following  notation  is  used.  Bold  lower-case  font 
indicates  a  vector,  while  bold  upper-case  font  indicates  a  matrix,  Pr(*)  indicates  the 

probability  of  an  event,  £{*}  indicates  the  expected  value,  (*)r  indicates  the  transpose 
operation,  (*)  indicates  the  Hermitian  transpose  operation,  ||*||  indicates  the  2-norm 
operation,  I  is  the  identity  matrix,  1/  indicates  a  column  vector  of  length  i,  tr(*) 
indicates  the  trace  of  a  matrix,  j  indicates  an  imaginary  quantity,  c  is  the  signal 
propagation  speed,  and  X  is  the  signal  wavelength. 

A.  WIRELESS  SENSOR  NETWORKS 

Advances  in  wireless  communications  coupled  with  low  cost  multifunctional 
sensors  have  given  rise  to  the  growth  of  WSNs.  A  WSN  is  a  network  comprised  of  many 
sensor  nodes  deployed  in  or  near  a  phenomenon  to  be  observed.  Typical 
implementations  have  identical  sensor  nodes  collecting  data  and  routing  it  through  a 
shared  channel  to  a  sink  node  for  processing/storage.  An  illustration  of  this  concept  is 
shown  in  Figure  1.  Often,  the  deployment  mechanisms  are  not  predetennined,  thus 
requiring  arbitrary  deployment  schemes  in  inaccessible  locations.  The  sensor  network 
nodes  collaborate  to  achieve  a  common  goal.  Some  examples  of  this  collaboration  can  be 
found  in  [25].  For  example,  sensor  nodes  can  be  distributed  across  a  geographic  region 
to  monitor  rainfall  to  support  flood  predictions,  covert  sensors  could  be  deployed  along 
known  supply  lines  to  monitor  enemy  troop  movements,  or  RFID  sensors  could  be 
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deployed  throughout  a  warehouse  for  inventory  control  [25].  A  few  other  notable 
research  areas  where  WSNs  have  been  applied  are  environmental  monitoring  [26], 
military  applications  such  as  acoustic  detection  of  helicopters  and  chemical  weapons 
detection  [27],  and  logistical  applications  such  as  global  shipment  tracking  [28]. 


1.  Ground-based,  Stationary  Wireless  Sensor  Network  System 
Architecture 

Sensor  networks  have  many  benefits,  but  the  realization  of  an  actual  network  can 
be  complex  with  many  design  factors  to  consider.  Some  key  factors  are  fault  tolerance, 
scalability,  production  cost,  hardware  constraints,  sensor  network  topology,  operating 
environment,  transmission  media,  power  consumption,  and  protocol  stack  [29]. 

A  major  aspect  of  WSNs  is  the  hardware  components  of  each  sensor  node.  A 
typical  sensor  node  has  four  basic  components:  a  sensor,  processor,  transceiver,  and 
power  source.  There  can  also  be  some  application  dependent  components,  such  as 
localization,  power,  and  mobility  systems.  Each  system  plays  a  key  role  in  this  body  of 
research.  The  key  components  of  such  a  sensor  node  are  shown  in  Figure  2. 


8 


Figure  2.  Key  components  of  a  sensor  node. 


2.  Network  Time  Synchronization 

Clock  synchronization  is  an  underling  requirement  for  many  WSN  applications. 
A  key  support  function,  time  synchronization  is  vital  to  such  applications  as  data  fusion, 
transmission  time  synchronization,  time-based  channel  scheduling,  and  node  sleep/wake 
scheduling  [30].  Most  studies  on  clock  synchronization  have  revolved  around  protocol 
design  [31],  but  at  its  core  clock  synchronization  is  a  parameter  estimation  problem,  thus 
meriting  a  signal  processing  approach.  The  clock  in  each  WSN  node  can  be  defined  as 


c(t)  =  t. 


(1) 


where  t  is  the  ideal  reference  time.  Realistic  systems  deviate  from  the  ideal  clock  due  to 
imperfect  clock  oscillators,  resulting  in  the  general  clock  function 


c{t)  =  sct  +  oc. 


(2) 


where  sc  is  the  clock  skew  and  oc  is  the  offset.  From  (2),  we  can  see  that  if  not 
addressed,  clock  error  can  grow  linearly  over  time.  For  example,  a  typical  crystal-quartz 
oscillator  commonly  used  in  WSNs  can  have  frequencies  that  vary  up  to  40  ppm,  causing 
potential  clock  errors  up  to  ±  40  /vs  per  second  [30].  Although,  clock  synchronization  is 
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a  concern  in  many  WSN  applications,  this  research  does  not  address  this  problem  and 
assumes  perfect  synchronization  across  the  network. 


3.  Elevated,  Mobile  Wireless  Sensor  Networks 

It  is  well  known  that  placing  a  received  antenna  on  an  elevated  platform  improves 
signal  reception.  By  elevating  the  antenna,  the  communication  link  is  less  likely  to  be 
obstructed  while  promoting  a  line-of-sight  signal  path.  From  [32],  the  received  power  at 
an  antenna  over  a  reflective  surface  can  be  expressed  as 


KK 

d2 


(3) 


where/]  is  the  transmitted  power,  hr  is  the  receiver  antenna  height,  ht  is  the  transmitter’s 
antenna  height,  d  is  the  distance,  Gr  is  the  gain  of  the  receive  antenna,  Gt  is  the  gain  of 
the  transmit  antenna,  and  L0  is  other  associated  losses.  From  this  expression,  we  can  see 
that  an  increase  in  receiver  or  transmitter  height  significantly  decreases  the  signal  path 
loss,  i.e.,  increases  the  ratio  of  P  toPr  .  For  example,  given  a  scenario  in  which 

h.  —  ht  —  2  m ,  Gr  =  Gt  —  L{)  =0  dB  and  d  =  5000  m  ,  the  resulting  path  loss  is  136  dB. 
By  increasing  both  hr  and  ht  to  100  m,  the  path  loss  is  decreased  to  68  dB. 

Mobility  in  sensor  networks  has  many  advantages  with  regard  to  source 
localization  and  collection.  In  this  research,  it  is  not  the  velocity  of  a  node  that  is 
important  but  its  ability  to  reposition  and  maintain  position.  For  example,  in 
collaborative  beamforming,  a  mobile  sensor  network  can  be  reconfigured  to  achieve  a 
desirable  array  factor  [2],  [33].  In  hyperbolic  localization,  a  sensor  network  can  be 
reconfigured  to  obtain  optimal  sensor  to  target  formations  for  increased  location 
estimation  accuracy  [17],  [34]. 

To  leverage  the  advantage  of  both  elevation  and  mobility,  the  proposed  scheme  is 
developed  in  the  context  of  an  elevated,  mobile  WSN;  elevated  to  promote  a  LOS  signal 
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path  for  increased  received  power  and  range,  and  mobile  for  control  over  sensor-target 
geometry  for  increased  localization  accuracy. 


4.  Multirotor  UAV  as  an  Elevated,  Nobile  WSN  Sensor  Node 

Sensor  nodes  for  WSNs  have  been  developed  for  different  operating 
environments,  such  as  stationary  ground  networks  [35],  ground  mobile  networks  [36], 
and  airborne  networks  [37].  For  signal  collection,  the  ideal  operating  environment  is 
elevated  above  physical  obstructions  and  noisy  surface  environments.  Such  a  vantage 
point  extends  communications  range  and  promotes  a  line-of-sight  signal  path.  For  this 
research,  we  propose  the  use  of  an  elevated,  mobile  WSN.  A  possible  realization  of  this 
network  is  via  the  use  of  multirotor  unmanned  aerial  vehicles  (UAVs)  [38],  [39].  These 
UAVs  can  provide  a  mobile  sensor  platform  that  in  contrast  to  fixed-wing  UAVs  can 
maintain  a  given  position.  This  research  does  not  undertake  sensor  network 
implementation  using  multirotor  UAVs;  however,  the  simulations  in  this  dissertation 
utilize  the  concept  of  an  elevated  and  mobile  sensor  node  realized  by  multirotor  UAVs. 

Although  the  modern  multirotor  UAV  has  only  recently  been  developed,  it 
remains  a  relatively  simple  machine.  For  instance,  a  quadrotor  has  four  rotors  held 
together  with  a  rigid  frame  [40].  Control  of  such  vehicles  is  basically  accomplished 
through  the  differential  manipulation  of  each  rotor’s  thrust  [38],  [39].  In  contrast  to 
fixed-wing  UAVs,  multirotor  UAVs  are  under-actuated,  with  the  remaining  degrees  of 
freedom  controlled  through  system  dynamics.  Having  demonstrated  their  effectiveness 
in  many  applications,  such  as  environment  mapping  and  monitoring  [41],  [42], 
transportation  and  construction  [43],  and  wireless  communication/networking  [44],  they 
have  become  a  popular  research  platform.  For  more  information  regarding  multirotor 
UAVs,  see  [38],  [39]. 

a.  Multirotor  UA  V  Station  Keeping 

There  are  many  aspects  to  multirotor  UAVs  to  consider  when  used  as  a  sensor 

node.  Of  key  concern  is  its  ability  to  station  keep  a  given  position,  for  any  deviation 

from  its  given  position  due  to  system  errors  and/or  wind  results  in  sensor  position  errors. 

The  station  keeping  performance  of  a  standard  proportional-integral-derivative  (PID) 
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controlled  quadrotor  in  the  presence  of  wind  is  found  in  Figure  3.  Here,  the  wind  is 
modeled  using  a  Dryden  model  [45].  From  the  illustration,  we  can  see  the  observed 
position  errors  due  to  station  keeping  operations  do  not  easily  conform  to  any  standard 
distribution.  The  actual  position  errors  of  a  multirotor  UAV  are  mainly  governed  by 
wind,  GPS  accuracy,  and  flight  controller  scheme  [46],  [47],  [48].  In  the  absence  of  such 
knowledge,  we  resort  to  two  basic  position  error  distributions,  uniform  and  Gaussian. 


a) 


!>) 


Figure  3.  Simulated  quadrotor  station  keeping  in  the  presence  of  Dryden  modeled 
wind,  a)  UAV  three-dimensional  (3D)  motion,  b)  histogram  of  two- 
dimensional  (2D)  position  errors,  from  [48]. 


B.  HYPERBOLIC  SOURCE  LOCALIZATION 

Hyperbolic  localization  is  the  process  of  determining  the  location  of  a  signal 
emitter  based  on  difference  of  a  signal’s  time-of-arrival  (TOA)  across  an  array  of  sensor 
nodes.  Because  of  its  dependence  on  the  difference  in  TOA,  it  is  also  called  localization 
via  time-difference-of-arrival  (TDOA)  measurements.  Hyperbolic  localization  has  been 
studied  extensively  and  can  be  found  in  many  applications  [49],  [50],  [51].  Hyperbolic 
location  is  commonly  applied  in  WSN  as  a  passive  technique;  passive  in  the  sense  that 
the  receiving  sensors  require  little  information  or  synchronization  with  the  transmitter. 
Two-dimensional  (2D)  hyperbolic  localization  begins  when  a  signal  is  received  by  three 
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or  more  synchronized  sensors  [5].  The  transmitter’s  location  is  then  determined  by  a  set 
of  hyperbolic  equations  that  are  derived  from  the  TDOA  measurements. 

To  derive  this  system  of  hyperbolic  equations,  let  us  consider  a  signal’s  TOA  at 
the  Ith  node  of  an  M  node  WSN  as 


ti=t0+di/c,  (4) 

where  to,  di,  and  c  are  the  signal’s  time-of-emission,  distance  from  the  transmitter  to  the 
7th  node,  and  signal  propagation  speed,  respectively.  Typically,  these  TOA  measurements 
are  obtained  using  the  cross-correlation  method  [16].  By  taking  the  difference  between  a 
pair  of  nodes’  measured  TOA,  we  eliminate  t0  and  can  obtain  a  set  of  TDOA  equations  as 

ti+l-tl=(di^-d])/c,  (5) 

with  /  =  2,3,...,  M  .  This  can  be  rewritten  to  generate  hyperbolic  equations  that  include 
the  transmitter’s  and  the  7th  sensor’s  location  in  Cartesian  coordinates  as  [52] 

c Om  - 0  =  V(x'+ 1 “ xt)  +  U-+i  - y, j  - v(*i - ■ x7 )  +  (ki - y, )  (6) 

where  ( xt , yt )  is  the  7th  node’s  location  and  ( xt,yt )  is  the  transmitter’s  location.  With  a 

minimum  of  three  sensor  nodes  yielding  two  simultaneous  equations  in  the  form  of  (6),  a 
two-dimensional  location  estimate  can  be  determined  [5]. 

Since  the  system  of  equations  governing  this  reconstruction  is  non-linear,  solving 
them  becomes  a  non-trivial  matter.  There  are  two  ways  to  deal  with  this  issue.  One  is  to 
avoid  the  non-linear  equations  altogether  by  selectively  setting  the  coordinate  system 
[34],  [5],  and  the  other  linearizes  the  function  through  Taylor  series  expansion  about  a 
reference  point  [17],  [52].  Additionally,  the  research  in  this  dissertation  is  focused  on 
elevated,  mobile  networks  with  an  emphasis  on  a  line-of-sight  (LOS)  signal  path  and 

does  not  address  the  non-LOS  case.  For  non-LOS  solutions,  see  [53]  and  [54], 
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1.  Weighted  Least-Squares  Estimation 

The  least-squares  method  is  the  standard  approach  for  parameter  estimation,  given 
a  set  of  measurements  that  are  derived  from  a  linear  function  of  inputs  and  parameters 
[55].  The  least-squares  solution  is  obtained  as  a  set  of  ^parameters  that  minimizes  the 
sum  of  the  squared  errors: 


ELS=^(yt-f(xnmLs))2  >  (7) 

1=1 

where  Nm  is  the  number  of  noisy  measurements,  /  (x,ft>/  s.  ]  is  a  linear  function  of  the  7  th 
input  variable  xj  and  the  Nm  x  1  solution  parameter  vector  coLS ,  and  yt  is  the  7th  noisy 
measurement  defined  as 


yi=f(xi,(oLS)+£w,  (8) 

with  sw  being  the  zero  mean  Gaussian  measurement  error.  It  is  important  to  clarify  the 
linear  nature  of  /(x;,toi5).  This  function  is  considered  a  linear  model,  when  the 

relationship  between  the  input  variable  and  the  output  is  a  linear  parameterization  of  (0LS. 
When  the  measurements  are  taken  at  different  levels  of  uncertainty,  the  weighted  least- 
squares  approach  has  been  shown  to  be  optimal  [56],  [57].  The  weighted  least-squares 
solution  minimizes  the  sum  of  weighted  squared-errors: 

Nm  2 

^wls  =  ^  wLs,i  (k,  ~  f  (xi 5 0Jm.s ))  >  (9) 

7=1 

where  wLS  .  is  the  7th  weight  and  (0WLS  is  the  N o)  x  1  parameter  vector.  The  weighted 
least-squares  estimate  in  matrix  form  is  expressed  as  [56] 
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awls=(xTWX)'1  XTWy,  (10) 

where  y  is  Nm  x  1  vector  of  noisy  measurements,  the  Nm  x  Nm  weighting  matrix  is 


where  cr2  is  the  measurement  variance  and  the  N (0  x  Nm  input  variable  matrix  is  [56] 

df{XVMLS)  _  Qf(XVt0Ls) 
dcot  dwN 

df{xNm,a>Ls)  df(xNm,(oLS) 

doix  dct)N 

2.  Maximum-Likelihood  Estimation 

The  maximum-likelihood  estimator  also  requires  the  system  of  equations  to  be 
linear.  The  term  “maximum-likelihood”  is  used  because  the  solution  maximizes  the 
likelihood,  i.e.,  the  statistical  model  of  the  estimate  matches  that  of  the  measurements. 
Since  the  maximum-likelihood  estimator  is  both  asymptotically  unbiased  and  efficient, 
i.e.,  achieves  the  Cramer-Rao  lower  bound  (CRLB)  [15],  [58],  it  has  become  widely 
adopted  in  the  field  of  parameter  estimation. 

To  derive  the  maximum-likelihood  solution,  we  consider  an  Mxl  noisy 
measurement  vector  given  by 

j  =  ^((o,jr)  +  £,  (13) 
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where  to  is  the  N(a  x  1  vector  of  unknown  but  nonrandom  set  of  parameters  to  be 
estimated,  z(co,x)  is  a  function  of  to  and  the  input  vector  x,  and  £  is  the  zero  mean 
Gaussian  measurement  error.  The  likelihood  ofj  for  a  given  to  is  governed  by  its 
conditional  probability  density  function  and  is  expressed  as  [52] 


f*.(y  l<»)= 


o*r¥. 


(  T  \ 

exp  -(1  /  2)  [j  -  z  (to,  x)]  C;1  [ y  -  z  ( a),  x)]  , 


(14) 


where  •  denotes  the  determinant  of  a  matrix,  and  C  is  the  covariance  matrix  of  the 
measurement  error  s  and  is  defined  as 


C  = 


(15) 


Finally,  the  maximum-likelihood  estimator  of  to  can  be  obtained  by  minimizing  the 
exponent  quadratic  of  (14)  [52] 

Qml  =  [. y  -z((o,x)J  c;1  [y  -z(co,  x)]  . 

We  now  encounter  the  same  problem  as  with  the  weighted  least-squares  estimator,  where 
z(cu,x)is  nonlinear,  i.e.,  the  relationship  between  y  and  .v  is  not  a  linear  parameterization 
of  to .  The  standard  solution  is  then  to  linearize  the  functions  through  a  Taylor  series 
expansion  about  a  reference  point  cor .  Using  only  the  first  two  terms  of  the  expansion, 
we  have  the  following  approximation  [52] 

zT  (oj,x)  «  z(to,x)  +Gts  ( (0-60 '),  (16) 
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where  Gts  is  the  M  x  NM  gradient  matrix  given  by 


Gts  = 


dz. 


dco, 


dz 


M 


dco, 


dzx 


dco v 


dz 


M 


dco. 


(17) 


The  maximum-likelihood  estimator  can  then  be  obtained  as  [52] 

<  =cor+(GlC-lGtsyGlC-l(y-z{cor,x)).  (18) 

It  is  important  to  note  that  the  maximum-likelihood  estimator  can  be  iteratively 
processed  if  there  are  multiple  sets  of  yTD  .  These  iterations  can  be  used  to  overcome  the 
effects  of  an  inaccurate  cor  or  compensate  for  increased  measurement  noise.  The 
iterative  performance  of  a  15  node  maximum-likelihood  estimator  initialized  with  an 
inaccurate  cor  is  shown  in  Figure  4.  With  co  set  1100  meters  away  from  the  true 
location,  we  can  see  that  the  estimator  quickly  approaches  the  asymptotic  optimal 
minimum  error  value  [59]. 


Inte  rations 


Figure  4.  Iterative  performance  of  a  hyperbolic  maximum-likelihood  estimator  in 

the  presence  of  an  inaccurate  cor  estimate. 
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3.  Optimal  Formations  for  Source  Localization 

To  derive  the  optimal  formations  for  hyperbolic  localization,  we  first  determine 
its  Cramer-Rao  lower  bound  (CRLB)  [17].  For  localization  using  polar  coordinates,  the 
CRLB  is  derived  using  the  Fisher  information  matrix  of  the  hyperbolic  localization.  The 
Fisher  infonnation  matrix  is  given  by  [17] 


0(<y) 


Yf  Zf 


dd°  dd° 

dr  d6 


dd°  del0 
dr  d6 


(19) 


where  r  is  range,  0\s  bearing,  d°  is  the  ( M  - 1  ]  x  1  true  TDOA  vector,  and  the 
(M-l)x(M-l)  covariance  matrix  of  the  TDOA  measurement  error  Ce  is  expressed  as 
[17] 


C  = 


(20) 


where  o\  is  the  range  measurement  error  variance,  I  is  the  ( M  —  1 )  x  ( M  - 1 )  identity 
matrix,  and  lM_j  is  a  column  vector  of  (M  —  l)  ones.  The  elements  of  the  Fisher 
information  matrix  from  (19)  are  [17] 


I,  =■ 


2  Mcr2Rc2rt4 


M 

M£r>n4(0, -e„)~ 


1=1 


M 


i= 1 


(21) 


Y  = _ - _ 

F  ,  f  J2  2  2 


M cr-c  rt  I  ,=1 


M 

M^r?sm3(0t-0.)- 


M 

£r2  sin2  (0,-3) 


1=1 


M 

£r.sin  [9,-9) 


i= 1 


(22) 


and 
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Male1 


sin2  (9,-9,)- 


i= 1 


^r.sin  (9,  -9,) 


Z=1 


(23) 


where  c  is  the  signal  propagation  speed  and  (r,  6()  is  the  /lh  sensor’s  polar  location. 

Based  on  the  use  of  the  maximum-likelihood  estimator,  which  is  asymptotically  efficient, 
the  objective  is  to  choose  a  sensor  formation  that  will  decouple  the  range  and  bearing 
estimates,  i.e.,  a  formation  that  makes  the  element  YF  =  0.  In  [17],  it  is  shown  that  a 
system  of  concentric  circles  satisfies  this  condition.  From  within  this  concentric  circle 
formation,  we  can  select  a  sub-fonnation  that  either  maximizes  XF  or  ZF,  i.e., 
minimizes  the  CRLB  of  range  or  bearing 


CRLB 


Mr 4 


O  2  l  4 

8 <rRc  rt 


-sin  9, 


and 


n  _J__ 

^CRLB  ry 


2Mra  •  2  . 

2~^rsm  9, 

<JRC 


(24) 


In  [52],  Torrieri  showed  that  the  fonnation  of  the  sensors  in  relation  to  the  signal 
emitter  affects  the  localization  accuracy.  This  degradation  in  accuracy  is  termed  the 
geometric  dilution  of  precision  (GDOP).  For  the  maximum-likelihood  estimator,  it  is 
given  by  [52] 


where 


ML 


=  tr 


(GlC~slGtsy  )/{cas), 


1=1 


(25) 


(26) 
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is  the  average  variance  of  the  arrival  times  with  cr,2(  being  the  variance  of  the  TOA 

measurements  at  the  7th  node.  To  combat  this  degradation,  we  employ  optimal  sensor 
formations  that  minimize  the  effects  of  GDOP  [17].  Through  inspection  of  the  CRLB,  it 
can  be  shown  that  the  sensor  formations  from  Figure  5  minimize  GDOP  [34],  [17]. 

Each  of  the  three  formations  in  Figure  5  is  suited  to  different  applications 
depending  on  the  estimated  parameter  of  interest.  Here,  J  is  the  number  of  sensor  groups 
and  M  is  the  total  number  of  sensors.  The  circular  formation(j  =  3  or,/ >5)  provides 
optimal  bearing  and  range  estimates,  independent  of  the  targets  bearing.  The  line 
formation  («/  =  2)  delivers  the  best  range  and  bearing  estimates  but  is  dependent  on  the 

target’s  bearing  [17],  [34],  The  cross  formation(  J  =  4)  delivers  performance  comparable 

to  the  circular  fonnation  with  minor  dependence  on  the  target’s  bearing.  Although  some 
of  the  fonnations  in  Figure  5  may  not  look  like  concentric  circles,  like  the  linear  or  cross 
formations,  they  are  indeed  two  concentric  circles.  One  large  circle  surrounding  a 
smaller  one  whose  radius  is  zero. 


Linear  (J  =  2) 


Cross  (J  =  4) 


Figure  5.  Optimum  array  geometries  for  position  estimation  (minimum 
uncertainty  area)  using  TDOA  localization.  J  =  number  of  sensor  groups,  M 
=  total  number  of  sensor,  from  [17]. 


For  more  information  regarding  these  optimal  formations  see  [17]  and  [34], 
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4.  The  Airborne  Symmetrical  Line  Array  Network  Configuration 

In  this  section,  we  describe  the  term  airborne  symmetric  line  array  (ASLA) 
formation.  The  ASLA  formation  consists  of  multiple  elevated,  mobile  sensors  in  a  line 
formation  symmetric  about  the  origin.  This  formation  is  based  on  the  linear  formation  in 
Figure  5  and  is  illustrated  in  Figure  6.  It  consists  of  three  sensor  stacks  that  contain 
Ns  =  M  /  3  sensors  each.  The  center  stack  is  located  at  the  origin,  with  two  outer  stacks 

flanking  the  center  at  a  distance  of  r  .  For  hyperbolic  localization,  this  formation  has 

been  show  to  minimize  the  maximum-likelihood  error  variance  [17].  Note  that  each 
sensor  in  a  stack  must  be  separated  by  1/2  to  guarantee  independent  noise  and  to 
minimize  mutual  coupling  [34], 


y-axis 


Figure  6.  Illustration  of  the  ASLA  formation. 


The  linear  formation  is  special  in  the  sense  that  it  delivers  the  best  performance  in 
bearing  and  range  estimates.  The  caveat  is  that  this  optimality  is  dependent  on  the  target 
emitter’s  bearing  [34],  A  Monte  Carlo  simulation  illustrating  the  dependence  on  the 
target  bearing  is  shown  in  Figure  7.  Here,  the  hyperbolic  localization  accuracy  of  an 
eight  node  ASLA  formation  is  compared  with  an  equal  numbered  circular  formation  at 
varying  target  bearing.  From  the  simulation,  we  see  that  the  bearing  of  the  target  has 
little  to  no  effect  on  the  circular  estimate,  whereas  the  linear  formation  performs  better 
when  the  target  bearing  is  approximately  ±  30  degrees  from  the  array’s  normal  direction. 
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Even  though  conditionally  optimal,  its  high  accuracy  is  the  key  reason  the  linear 
formation  was  chosen  for  the  proposed  scheme. 


Figure  7.  Accuracy  of  position  estimate  for  fixed  linear  formation  (blue)  and  fixed 

circular  formation  (red). 

5.  TDOA  Outlier  Measurement  Diagnostics 

With  the  proposed  localization  technique  implemented  on  an  elevated,  mobile 
sensor  network,  there  will  be  increased  measurement  noise  due  to  sensor  position  and 
other  sensor  network  related  errors.  This  increase  in  measurement  noise  creates  a  need 
for  methods  to  mitigate  the  effects  of  statistical  outliers  or  errors  on  parameter  estimation. 
In  this  section,  we  cover  some  statistical  tools  used  in  our  proposed  localization  scheme 
for  robustness  to  such  errors. 

To  assess  the  influence  of  the  7th  measurement/observation  in  a  least-squares 
solution,  we  can  calculate  the  solution  with  and  without  it.  Calculating  the  amount  of 

divergence  with  and  without  the  measurement  is  termed  “single  case  diagnostics”  [18]. 
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To  perform  this  diagnostic,  it  is  important  to  consider  which  measure  of  divergence  is 
appropriate  for  the  solution.  A  standard  approach  is  to  use  the  squared-Mahalanobis 
distance  given  by  [18],  [19], 

Dm  =(v,-£{v,})C„-'(v,-£{v,})r,  (27) 

where  the  vwis  the  Mx 2  matrix  containing  M  location  estimates,  v;  is  the  1x2  vector 

containing  the  z'th  location  estimate,  and  CJs  the  covariance  matrix  of  vm .  The  purpose 

of  the  Mahalanobis  distance  is  to  discern  how  much  influence  the  zth  measurement  has  on 
the  solution  given  the  covariance  of  each  of  the  measurements. 

To  illustrate  this  concept,  we  consider  a  set  of  200  weighted  least-squares  location 
estimates  shown  in  Figure  8.  In  this  set,  there  is  one  outlier  shown  in  green,  which  we 
wish  to  reject.  Using  the  Euclidean  distance  from  the  mean  as  a  measure  of  divergence, 
we  see  that  the  outlier  has  the  same  divergence  as  the  standard  estimate  shown  in  red.  In 
contrast,  using  the  Mahalanobis  distance  as  the  measure  of  divergence,  we  find  that  the 
outlier  has  a  divergence  of  34.5  where  as  the  standard  estimate  has  a  divergence  of  0.9; 
hence,  the  outlier  can  be  easily  identified  and  discarded. 
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Figure  8.  Scatter  plot  of  200  hyperbolic  location  estimates. 
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c. 


BEAMFORMING 


Beamforming  is  an  efficient  signal  processing  technique  that  enables  space- 
division  multiple  access  [6].  The  benefits  and  various  applications  of  traditional 
beamfonning  have  been  well  documented  [60],  [61].  The  extension  of  beamforming  to 
wireless  sensor  networks  has  given  rise  to  the  concept  of  collaborative  beamforming  [2], 
[62]  in  which  a  network  of  synchronized  and  distributed  sensors  is  used  to  fonn  a 
distributed  beamforming  array. 

Beamforming  begins  with  a  plane  wave  signal  impinging  on  the  array  from  a 
direction-of-arrival  9t.  With  the  signal  arriving  at  each  node  with  a  corresponding  phase 
delay,  beamforming  uses  phase  shifters  to  nullify  this  delay  and  constructively  sums  each 
node’s  signal.  The  concept  of  a  plane  wave  falling  across  a  uniform  linear  array  of  M 
nodes  where  the  phase  delay  is  a  function  of  the  distance  between  each  node  d  is 
illustrated  in  Figure  9. 


Array  normal 


Array  node  i  -  1  2  3  .  M 

Figure  9.  Plane  wave  impinging  on  a  uniform  linear  array. 
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1.  Phase  Shift  Beamforming 

With  the  phase  of  the  incoming  signal  at  node  i  =  1  set  to  zero,  the  combined 
signal  of  a  uniform  linear  array  with  M  nodes,  commonly  termed  the  array  factor  can  be 
expressed  as  [63] 


M 


4*(<?.0„)  =  £ 


we 


j(i-l)  /fdsinfy 


i= 1 


(28) 


where 


wi=Aiem™e-  (29) 

is  the  complex  weight  with  magnitude  At  and  phase  -/3dsmOt,  A  is  the  signal 
wavelength  and  /3  -  ( 2/rj  /  A .  From  (28),  it  is  apparent  that  the  array  yields  a  maximum 
response  in  the  direction  of  the  beam  steering  angle  0sa .  By  controlling  the  response  of 
the  array  as  a  function  of  0sa,  the  spatial  diversity  required  for  space-division  multiple 
access  is  achieved.  Along  with  this  diversity,  beamforming  also  provides  a  signal  gain 
equal  to  the  number  of  receiving  elements/nodes. 

There  are  many  parameters  that  dictate  the  resulting  array  factor  of  a 
beamforming  array.  Of  key  importance  to  this  research  is  the  number  elements/nodes 
and  their  placement  relative  to  each  other.  The  effect  of  increasing  the  number  of 
elements/nodes  in  a  uniform  linear  array’s  array  factor  is  illustrated  in  Figure  10.  The 
first  effect  is  the  narrowing  of  the  main  beam’s  beamwidth.  The  main  beam  is  defined  as 
the  lobe  that  contains  the  direction  of  the  array  factor’s  maximum  response.  The 
beamwidth,  also  called  the  half-power  beam  width,  is  typically  defined  as  the  angle  range 
at  which  the  main  beam  falls  above  3  dB  of  the  maximum  response.  The  relationship 
between  the  beamwidth  and  the  number  of  nodes  in  a  linear  array  is  given  by  [63] 
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(30) 


~  0.88/1 
w~  dM  ' 

The  second  effect  is  the  location  and  the  number  of  sidelobes  present  in  the  array  factor. 
Sidelobes  are  defined  as  any  lobe  other  than  the  main  beam  lobe.  For  uniform  linear 
arrays,  the  magnitude  of  the  nearest  sidelobe  to  the  main  lobe  is  on  average  13  dB  lower 
[63]. 


0 


Normalized  Power  (dB) 


Figure  10.  Array  factor  of  a  15  isotropic  node  linear  array  with  d  =  / i/2  . 

2.  Adaptive  Beamforming 

Adaptive  beamforming  is  the  process  of  adapting  an  antenna’s  array  factor  in  a 
manner  that  both  increase  the  gain  of  a  target  signal  while  decreasing  the  gain  of  all 
interfering  signals.  This  is  accomplished  by  placing  a  complex  weight  on  each  antenna 
signal.  These  weights  modify  the  phase  and  allow  for  the  coherent  amplification 
(reduction)  of  a  target  signal  (undesired  signals).  In  this  section,  we  will  illustrate  the 
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basic  phase  shift  beamfonner  [64],  [65]  and  the  adaptive  minimum  variance  distortionless 
response  (MVDR)  beamfonner  [66],  [67]. 

The  MVDR  beamfonner  combats  the  effects  of  interfering  signals  by  placing 
nulls  in  their  direction.  This  is  achieved  by  selecting  the  array  weights  that  minimize  the 
output  noise  variance  [63].  The  anay  output  is  expressed  as 

y(t)  -  wHs  +  wHu ,  (31) 

where  w  is  an  Mxl  vector  of  the  complex  weights  to  be  detennined,  u  is  an  Mxl 
vector  of  the  sum  of  all  interfering  signal  vectors  i)i  =[\,e,/3dsm0i 1)</sin6!]r 
represented  as 

u  =  T^hui(tH’  t-32) 

and  s  is  an  Mxl  the  signal  vector  represented  by 

s  =  s(t)v,  (33) 

with  s(t)  being  the  desired  signal  and  the  steering  vector  v  =  [l,el/3dsm0'  l)dsm0,  ]T . 

By  minimizing  the  noise  contributions  of  the  interfering  signals  in  the  output,  we  can 
determine  the  optimal  weights  wo .  The  output  variance  can  be  shown  as  [63] 

o2v  =wH Rsw  +  wH Ruw .  (34) 

Assuming  that  Rs  =  E  \ss!/ 1  and  R  =  E  |  uu  "  |  are  known,  we  see  that  the  minimization 

of  the  output  variance  is  reduced  to  the  minimization  of  wH Ruw ,  which  results  in  an 
expression  for  the  optimal  weights  [63] 


27 
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(35) 


To  illustrate  both  the  phase  shift  and  MVDR  beamfonning  concepts,  we  consider 
a  scenario  where  a  target  signal  with  two  interfering  signals  is  received  by  the  15  node 
uniform  linear  array  shown  in  Figure  10.  The  target  signal-of-interest  is  shown  in 
Figure  11,  with  a  direction-of-arrival  (DOA)  of  45  deg  and  a  signal  amplitude  of  one. 
The  two  interfering  signals  are  arriving  at  a  DOA  of  30  and  50  deg,  each  with  an 
amplitude  of  10. 


Figure  11.  Transmitted  signal-of-interest. 

The  phase  shift  beamformer’s  response  to  the  interfering  signals  and  target  signal 
is  shown  in  Figure  12.  We  can  see  the  two  interfering  signals  have  fully  corrupted  the 
target  signal.  In  contrast  to  the  phase  shift  beamformer,  the  adaptive  MVDR  beamformer 
response  is  shown  in  Figure  13.  From  this,  we  can  clearly  see  that  the  MVDR 
beamfonner  is  able  to  nullify  interference  and  successfully  recover  the  target  signal. 
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Figure  12.  Response  of  15  node  phase  shift  beamformer  with  two  10  dB  interfering 

signals  at  DOA  of  30  and  50  degrees. 


Time  (seconds) 

Figure  13.  Output  signal  of  the  MVDR  beamfonner  in  the  presence  of  a  target 
signal  at  0t  =  45  deg  and  two  interfering  signals  at  Gi  =  30  and  50  deg  . 

A  comparison  between  the  array  factor  of  the  MVDR  and  that  of  the  phase  shift 
beamformer  is  provided  in  Figure  14.  From  the  illustration,  we  can  see  the  MVDR’s 
array  factor  contains  two  strong  nulls  in  the  directions  of  the  interference,  whereas  the 
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phase  shift  beamformer  has  a  response  of  approximately -13  dB  and  -3  dB  in  the  same 
directions. 


Figure  14.  Comparison  between  the  MVDR  and  phase  shift  beamformers  in  the 
presence  of  a  target  signal  at  6t  =  45  deg  and  two  interfering  signals  at 

=  30  and  50  deg . 


3.  Error  Effects  in  Beamforming 

Two  key  errors  to  account  for  in  beamforming  are  pointing  errors  and  unintended 
sidelobe  response.  Pointing  errors  are  defined  as  the  difference  between  the  array’s 
actual  steering  direction  and  that  of  the  maximum  response.  The  standard  deviation  of 
the  pointing  error  is  given  by  [4] 


PE~  4m 


(36) 


where  cr  is  the  standard  deviation  of  the  phase  errors.  With  the  beamwidth  for  a 
uniform  linear  array  as  defined  in  (30),  the  pointing  errors  can  be  expressed  as 
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0.88/lcr, 


(37) 
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The  expression  for  the  expected  side  lobe  level  is  given  by  [63] 


A 


sl 


a(l-^)+(TAa 

Ma1  A2 


(38) 


where  a  is  the  error  free  main  beam  gain,  cr"a  is  the  variance  of  amplitude  error  and 


A2  = 


E[e‘ 


.  If  the  phase  errors  are  zero  mean  and  Gaussian,  then  A  -e  a&* . 


4.  Sidelobe  Control  via  Array  Tapering 

For  a  uniform  linear  array,  the  sidelobe  response  can  be  controlled  by  adjusting 
the  gain  of  each  node’s  response  [68].  This  process  is  tenned  tapering  and  is  akin  to 
windowing  in  digital  fdtering  [69],  [70].  Tapering  an  array  provides  a  tradeoff  between 
the  beamwidth  and  the  sidelobe  levels.  With  this  in  mind,  we  will  briefly  examine  the 
characteristics  of  four  tapering  methods:  Unifonn  [71],  Dolph-Chebyshev  [72],  Taylor- 
Kaiser  [73],  and  Binomial  [74]. 

a.  Tapering  Methods  for  Uniform  Linear  Arrays 

The  uniform  tapered  array  is  actually  a  non-tapered  array  as  all  nodes  are 
uniformly  weighted  by  one,  which  results  in  the  standard  phase  shift  beamformer  with  an 
array  factor  of 


M 

afl  ( e ) = fy(m-i)/Wsin* 

m= 1 


(39) 


This  tapering  method  makes  no  consideration  to  sidelobe  response  and  has  the  highest 
sidelobe  level  when  compared  to  other  methods. 
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The  binomial  tapered  array  is  unique  because  it  has  zero  sidelobes,  but  this  is  at 
the  cost  of  a  wider  beamwidth.  The  weights  of  an  M-node  binomial  uniform  linear  array 
are  calculated  from  the  binomial  coefficients  as  follows  [70] 


(/«)  = 


(M- 1)! 


m !  ( M  - 1  -  m ) ! 


m  =  0,1  -1 . 


(40) 


The  corresponding  array  factor  is  given  by  [70] 


r  1 


2cos 


V  J 


(41) 


From  this  expression,  we  can  see  that  the  array  factor  decreases  to  zero  as  6  approaches 
+  tc  . 


The  Dolph-Chebyshev  [72]  taper  method  is  best  applied  when  the  sidelobe 
response  at  all  angles  must  be  kept  below  a  specific  value  R/x  .  The  taper  weights  for  the 

Dolph-Chebyshev  uniform  linear  array  are  calculated  using  the  Dolph-Chebyshev 
window  transform  [70],  [75] 


Wcheb  {®k  ) 
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/  •  \ 
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M  cos  1 

R0  cos 
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cosh(M  cosh  1  (i?0)) 


i  -  0,1, ...,M  -1. , 


(42) 


where  the  scaling  factor  R(}  is  defined  as  [75] 


R0  = cosh 


cosh^lO^20)^ 


M 


(43) 


The  array  weights  wDC(m)  are  then  computed  as  the  inverse  discrete  Fourier  transform  of 
WCheh  ( 0)k ) .  The  corresponding  array  factor  for  an  array  with  an  even  number  of  nodes  is 
expressed  as  [70] 
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(44) 
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The  Taylor-Kaiser  tapering  method  is  similar  to  the  Dolph-Chebyshev  method  but 
has  an  exponentially  decreasing  sidelobe  response.  The  array  weights  can  be  computed 
as  [75] 


wm  (m)  =  /0  naTK \Jl-m2  / M 2 


(45) 


where  I0  is  the  zeroth-order  modified  Bessel  function  of  the  first  kind  and  aTK  is  the 

parameter  used  to  control  the  sidelobe  level  response.  The  Taylor-Kaiser  array  factor  can 
be  expressed  as  [70] 


AFT  {6) 
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sinh  Ja^ft2  -(MO / 2) 
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(46) 


b.  Comparison  of  the  Taper  Methods 

A  comparison  of  the  four  taper  methods  is  illustrated  in  Figure  15.  From  the 
results,  we  see  that  the  beamwidth  is  inversely  proportional  to  the  sidelobe  response.  We 
can  see  that  the  unifonn  tapering  achieves  the  narrowest  beamwidth  with  the  highest 
sidelobe  response,  whereas  the  binomial  tapering  has  the  widest  beamwidth  with  zero 
sidelobe  response.  The  Dolph-Chebyshev  has  the  unique  ability  to  maintain  a  given 
maximum  sidelobe  response,  while  the  Taylor-Kaiser  taper  show  a  better  sidelobe 
response  at  the  cost  of  a  slightly  wider  beamwidth. 
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Figure  15.  Beam  width  and  sidelobe  level  tradeoff  for  four  different  tapering 

methods  using  a  15  node  ULA.  The  beamwidth  “3  dB  |jne  js  shown  in  red, 

from  [70]. 


5.  Grating  Lobe  Control  via  Virtual  Filling 

In  collaborative  beamforming  where  the  nodes  are  separated  by  a  distance  greater 
than  XU,  grating  lobes  in  the  array  factor  will  be  a  significant  issue  [71].  Grating  lobes 
are  lobes  that  have  the  same  magnitude  as  the  main  lobe.  To  minimize  the  effects  of 
grating  lobes,  the  array  can  be  virtually  filled  [76].  The  virtual  filling  technique 
minimizes  the  grating  lobes  by  filling  in  any  internode  distance  dv  greater  than  XU 
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with  virtual  nodes.  These  virtual  nodes  consist  of  interpolated  complex  signal  data 
derived  from  the  real  nodes.  This  filling  is  done  in  a  manner  that  turns  the  array  into  a 
virtually  filled  uniform  linear  array,  which  is  inherently  free  of  grating  lobes  [76]. 

A  single  snapshot  of  the  complex  signal  received  at  the  nth  node  in  a  collaborative 
beamforming  network  can  be  expressed  as 

nr 

A  =5X,-  exp ( j(/)ti ) exp ( -jfixn  sin6> ,.),  n  =  1,...,M  ,  (47) 

1=1 

where  N R  is  the  number  of  signals  impinging  on  the  array,  xn  is  the  «th  node’s  location 
on  the  x-axis,  Vti,  (j)ti,  and  6ti  are  the  z'th  signal’s  magnitude,  phase,  and  bearing, 
respectively.  This  expression  can  be  written  in  matrix  fonn  as  [76] 

EnVy=An,  (48) 

where  An  is  an  Mxl  vector  containing  the  complex  signal  snapshot  from  each  node,  En 
is  an  Mx  NR  matrix  given  by 

exp  ( -jpx,  sin  0tX )  •  •  •  exp  ( -j fix x  sin  0i  Nr  ) 

E„=  (49) 

exp ( sin 0t  l )  •••  exp(-y^xw  sin^  v,; ) 

and  the  NR  x  1  vector 

K,i  exP(M,i) 

vv=  :  .  (50) 

K.n,  exp(.M,  vs) 

Assuming  that  each  signal’s  bearing  is  known,  an  estimate  of  their  magnitude  and  phase 
can  be  obtained  using  the  least-squares  approach  expressed  as  [77] 
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(51) 


Vr=(E:E„fElA„. 

These  signal  estimates  are  then  used  to  represent  each  virtual  nodes  signal  output.  The 
output  of  each  virtual  node  is  expressed  as 

nr 

A  =  cxpf.M.i)cxP(^T  sin 0ti ),  v  =  1  ,  (52) 

i=l 

where  xv  is  the  position  of  the  vth  virtual  node  along  the  x-axis  and  Mv  is  the  number  of 
virtual  nodes  equal  to  dv  /  (A  /  2) . 

In  this  chapter,  we  covered  hyperbolic  localization  and  beamforming  within  the 
context  of  wireless  sensor  networks.  With  regard  to  localization,  we  presented  two 
localization  techniques  and  their  application  to  an  optimal  ASLA  sensor  formation.  We 
also  covered  the  use  of  outlier  detection  techniques  to  mitigate  measurement  noise.  The 
discussion  regarding  beamforming  focused  on  techniques  to  increase  array  robustness 
against  various  steering  errors  and  sidelobe  interference.  Together,  these  topics  provide 
the  necessary  details  to  present  our  proposed  signal  collection  scheme  in  the  next  chapter. 
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III.  SOLUTION  APPROACH 


This  chapter  begins  with  an  explanation  of  the  operating  scenario  and 
assumptions  used  in  this  dissertation.  We  provide  an  overview  of  the  proposed  scheme 
for  radio  frequency  source  localization  and  signal  collection.  We  then  define  the  figures 
of  merit  used  to  evaluate  the  scheme’s  perfonnance. 

A.  PROPOSED  SCHEME 

The  intended  operational  scenario  for  the  proposed  scheme  is  shown  in  Figure  16. 
With  a  distant,  stationary,  sporadic  or  repeating  source  emitter  located  at  (x,,y(),  there  is 

a  network  of  sensors  deployed  with  the  goal  of  locating  the  emitter  and  collecting  its 
signals.  In  this  scenario  and  in  this  dissertation  research  in  general,  we  assume  the 
following.  Sensor  node  deployment  is  restricted  to  a  circular  area,  each  node  has  an 
isotropic  antenna  with  matched  polarization,  and  all  nodes  are  synchronized  to  a  common 
clock. 

The  network  of  M  sensor  nodes  is  initially  deployed  into  an  ASLA  formation 
(introduced  in  Chapter  II)  where  they  form  an  ad  hoc  wireless  sensor  network.  Each 
node  consists  of  a  multirotor  UAV,  used  to  realize  the  concept  of  an  elevated,  mobile 
sensor  network.  The  ASLA  formation  shown  in  Figure  6  places  the  M  sensor  nodes  into 
three  node  stacks,  each  containing  an  equal  number  of  nodes.  These  three  stacks  are 
deployed  along  the  x-axis,  one  at  the  origin  and  the  other  two  flanking  the  center  on  both 
sides  at  a  distance  of  ra .  For  source  localization,  this  formation  is  used  to  minimize  the 
estimate  error  variance  [17],  [34],  For  signal  collection,  the  node  stacks  allow  for  robust 
signal  estimation. 
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Figure  16.  Operational  concept  of  an  elevated,  mobile  sensor  network  deployed  in 

an  ASLA  formation. 

A  schematic  diagram  of  the  proposed  source  localization  and  signal  collection 
scheme  is  shown  in  Figure  17.  The  scheme  begins  in  the  localization  phase  once  a  signal 
has  been  detected.  This  generates  a  set  of  TDOA  measurements  that  is  fed  to  the  location 
estimator.  Using  this  initial  location  estimate,  we  reorient  the  sensor  network  to  be 
perpendicular  to  the  target  emitter.  The  objective  of  this  network  orientation  is  to 
maximize  the  accuracy  of  the  location  estimate,  which  we  call  the  refined  estimate.  With 
the  refined  estimate  obtaining  an  improved  location  estimate,  this  information  is  used  in 
the  signal  collection  phase.  In  this  phase,  each  node  in  the  network  samples  the  signal 
and  transmits  these  samples  to  the  network’s  sink  node.  Located  in  the  center  of  the 
network,  the  sink  node  uses  a  combination  of  beamform  processing  and  signal  estimation 
to  combine  and  amplify  the  signal  samples  coherently.  This  makes  the  signal  collection 
phase  essentially  a  collaborative  beamforming  effort. 
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Figure  17.  Process  flow  for  propose  signal  collection  scheme. 


1.  Source  Localization 

The  objective  of  this  phase  is  to  obtain  a  precise  location  estimate  of  the  target 
signal  emitter.  To  accomplish  this,  we  employ  the  use  of  an  unbiased  and  efficient 
maximum-likelihood  location  estimator  [52].  Although  an  effective  location  estimator,  it 
requires  an  initial  position  estimate,  in  order  to  avoid  local  minima  solutions  [59].  To 
satisfy  this  requirement,  we  precede  the  maximum-likelihood  estimator  with  another 
estimator  that  does  not  require  initialization.  For  this  task,  we  propose  a  weighted  least- 
squares  estimator. 

Furthermore,  after  obtaining  the  initial  estimate  we  take  advantage  of  the 
network’s  mobility  and  re-orient  the  network  nodes  to  create  an  optimal  sensor  to  target 
geometry  [17].  This  new  orientation  minimizes  the  geometric  dilution  of  precision  [52], 
thus  minimizing  the  maximum-likelihood  estimator’s  error  variance.  Because  of  this 
further  refinement  of  the  location  estimate,  we  term  the  maximum-likelihood  estimator 
the  refining  estimator.  This  two-stage  localization  technique  is  illustrated  in  Figure  18. 
Starting  in  the  initial  configuration  shown  in  Figure  18  (a),  we  obtain  a  location  estimate 

( x, ,  v, )  which  is  used  to  calculate  the  emitter  bearing  estimate  0t.  Using  0! ,  we  reorient 
the  network  nodes  perpendicular  to  the  target  signal  emitter,  resulting  in  the  reoriented 


network  configuration  as  seen  in  Figure  18  (b). 
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Figure  18.  Network  reorientation  operation  for  the  two-stage  localization  technique. 


In  conjunction  with  both  estimators,  we  propose  to  implement  a  robust 
measurement  outlier  rejection  process  to  detect  and  reject  erroneous  TDOA 
measurements.  This  process  is  used  to  increase  the  localization  phase’s  robustness  to 
measurement  error  and  sensor  position  errors. 

2.  Sensor  Node  Position  Error 

In  this  research,  the  position  information  of  each  node  is  assumed  known. 
Fluctuations  in  each  node’s  position  due  to  their  station  keeping  operations  are  also 
assumed  [47],  [48],  [78].  The  effects  of  these  position  errors  are  the  main  concern  in  this 
dissertation.  With  each  node  in  the  network  envisioned  as  a  multirotor  UAV,  the  actual 
position  errors  are  mainly  governed  by  wind,  GPS  accuracy,  and  flight  controller  scheme 
[46]-[48].  In  the  absence  of  such  knowledge  we  resort  to  two  basic  position  error 
distributions,  uniform  and  Gaussian.  More  specifically,  we  model  these  small 
fluctuations  in  the  x  and  y  positions  as  two  independent  and  identically  distributed 
uniform  or  Gaussian  random  variables  8X  and  8  ,  respectively. 

For  the  uniformly  distributed  position  error  case,  the  probability  density  functions 
are  given  by 
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where  8  can  be  interpreted  as  the  maximum  position  error  in  the  x  or  y  direction  [48]. 

For  the  case  of  Gaussian  position  error,  both  random  variables  are  zero  mean  with 

2 

a  variance  of  cr .  The  probability  density  functions  are  expressed  as 
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where  cr  is  the  position  error  variance  in  the  x  or  y  direction  [48]. 

3.  Collaborative  Signal  Collection 

After  the  source  localization,  the  objective  of  the  signal  collection  is  to  collect  and 
amplify  the  target  signal.  The  main  concern  is  that,  while  the  mobility  of  each  node  is 
allowed  for  optimal  formations  in  the  localization  phase,  it  will  cause  array  phase  errors. 
To  achieve  the  objective,  despite  the  presence  of  such  errors,  we  propose  a  combination 
of  collaborative  beamforming  and  signal  estimation. 

In  this  approach,  we  exploit  the  sensor  grouping  of  the  ASLA  formation.  The 
ASLA  formation  contains  three  groups  of  sensors  placed  symmetrically  along  the  x-axis 
(see  Figure  6).  With  each  node  in  a  group  providing  a  noisy  sample  of  the  same  signal,  a 
sample  mean  can  be  calculated.  Using  each  group’s  sample  mean  in  a  collaborative 
phase  shift  beamformer,  we  reconstruct  and  amplify  the  signal. 

Furthermore,  since  the  ASLA  formation’s  array  factor  contains  grating  lobes  due 

to  its  large  inter-stack  distance  r,  it  is  extra  sensitive  to  interfering  signals.  To  increase 

the  robustness  of  signal  collection  against  such  signals,  we  employ  the  use  of  virtual 
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filling  [76].  This  array  processing  technique  is  used  to  manage  the  array’s  inherent  side 
lobes  and  grating  lobes  in  order  to  minimize  the  effects  of  interfering  signals. 

It  is  important  to  note  that  while  our  scheme  is  not  a  traditional  adaptive 
beamformer,  it  accomplishes  the  same  task.  Adaptive  beamforming  is  the  process  of 
manipulating  an  array’s  weights  in  order  to  amplify  a  target  signal  while  suppressing 
undesirable  signals.  Although  our  approach  uses  fixed  uniform  weights,  it  amplifies  a 
target  signal  using  array  rotation  and  suppresses  undesirable  signal  through  the  use  of 
virtual  filling  and  tapering  techniques. 

B.  PERFORMANCE  METRICS 

With  this  dissertation  heavily  focused  on  parameter  estimation,  it  is  critical  to 
define  perfonnance  metrics  as  a  means  to  support  and  analyze  the  proposed  theory  and 
scheme.  The  perfonnance  metrics  used  throughout  this  dissertation  is  defined  as  follows. 

For  localization  schemes,  a  widely  used  performance  metric  is  the  Cramer-Rao 
lower  bound  (CRLB)  [58].  The  CRLB  is  the  theoretical  minimum  solution  variance  an 
estimator  can  achieve.  Any  unbiased  estimator  that  can  achieve  the  CRLB  is  said  to  be 
efficient  [59].  For  our  proposed  network,  the  CRLB  for  location  (minimum  location 
uncertainty)  and  bearing  (see  Chapter  II)  are  expressed  as 

(55) 
and 

(56) 

where  M  is  the  number  of  sensor,  ra  is  the  distance  between  the  outer  and  center  node 
stacks,  <j2r  is  the  ranging  measurement  error  variance,  and  6t  is  the  target  source  bearing. 

For  estimation  processes,  the  root  mean-square  error  is  a  standard  perfonnance 
metric.  For  source  localization,  we  focus  on  two  root  mean-square  enor  values,  the  first 
being  root  mean-square  enor  of  location  estimates  defined  as 
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(57) 
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loc 


where  cot  is  a  2x1  vector  containing  the  true  source’s  location,  d)t  is  the  estimate  cot ,  and 
nt  is  the  number  of  trials.  The  second  value  is  the  root  mean-square  error  of  bearing 
estimates  defined  as 


4  = 


(58) 


where  0t  is  the  true  source  bearing  and  0t  is  its  estimate. 

For  signal  estimation,  we  focus  on  the  root  mean-square  error  of  magnitude  and 
phase.  The  root  mean-square  error  of  magnitude  estimates  is  defined  as 
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(59) 


where  Vt  is  the  true  signal  magnitude  and  Vt  is  its  estimate.  The  root  mean-square  error  of 
phase  estimates  is  defined  as 


4 


(60) 


where  </>t  is  the  true  signal  magnitude  and  (f>t  is  its  estimate. 

To  quantify  the  error  in  a  normalized  array  factor  with  and  without  position 
errors,  we  define  the  normalized  array  factor  error  as 


(61) 
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where  Ot  is  the  signal  bearing,  AF  is  the  normalized  error  free  array  factor,  and  Ae  is  the 
nonnalized  array  factor  with  position  errors. 

In  this  chapter,  we  provided  an  overview  of  the  proposed  scheme  for  source 
localization  and  signal  collection.  We  described  the  use  of  two  sequential  location 
estimates  that  support  a  follow  on  signal  collection.  We  then  described  a  set  of 
performance  metrics  used  to  analyze  the  scheme’s  perfonnance.  In  the  next  chapter,  we 
present  a  two-stage  localization  technique  capable  of  approaching  the  CRLB  and  propose 
a  measurement  outlier  rejection  process  to  increase  localization  robustness. 
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IV.  ROBUST  TWO-STAGE  SOURCE  LOCALIZATION  FROM  AN 
AIRBORNE  SYMMETRIC  LINE  ARRAY  NETWORK 


In  this  chapter,  we  introduce  a  two-stage  source  localization  technique.  We 
examine  its  perfonnance  in  the  presence  of  measurement  and  sensor  position  errors.  To 
combat  these  effects,  we  propose  a  measurement  outlier  rejection  technique.  This 
technique  is  used  to  identify  and  reject  specious  TDOA  measurements  in  an  effort  to 
increase  the  localization  robustness.  We  then  use  simulations  to  support  the  theoretical 
development  and  analyze  the  perfonnance. 

A.  TWO-STAGE  HYPERBOLIC  LOCALIZATION 

The  objective  is  to  estimate  the  location  of  a  single  source  emitter  located  at 
(xt,yt).  The  sensor  network  is  deployed  in  the  initial  ASLA  configuration,  as  seen  in 

Figure  18  (a).  In  this  configuration,  the  network  detects  an  incoming  emission  and 
generates  a  set  of  TDOA  measurements.  Using  these  measurements  in  the  weighted 
least-squares  estimator  introduced  in  Chapter  II,  we  obtain  a  location  estimate.  This 
initial  localization  is  a  closed-form  estimator  that  does  not  require  any  initialization  [5], 
[79]. 


1.  Hyperbolic  Localization  via  Weighted  Least-Squares 

The  estimator  described  in  this  section  is  derived  from  the  linear  array  estimators 
found  in  [5]  and  [79],  and  a  generic  linear  array  formation  is  adapted  for  use  with  an 
ASLA  formation. 

Since  weighted  least-squares  estimation  is  primarily  for  linear  models,  its 
applicability  to  hyperbolic  estimation  is  not  immediately  obvious.  To  apply  the  weighted 
least-squares  approach  to  localization,  we  must  rewrite  (6)  to  obtain  a  set  of  linear 
equations.  To  do  this,  we  define  the  noisy  range  difference-of-arrival  measurement  of  a 
sensor  pair  as  [5],  [79] 


dij,  ~  d \k  +  si,k  ’ 


i,k  =  , 


(62) 
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where  ct-  k  is  the  true  range  difference  between  sensors  i  and  k  and  si  k  is  the  zero  mean 

Gaussian  range  difference  measurement  error.  Then  the  range  from  the  /lh  sensor  to  the 
source  emitter  is  determined  as  [5] 


rt  =yl(xi-xl )2+  (yt- ytf  .  (63) 

By  squaring  both  sides  of  (63),  we  obtain 

r;  =  Kt  - 2xtx,  - 2y,y ,+xf+yf,  i  =  (64) 

where  K.  =  xf  +  yf  and  r.  x  =  cci ,  =  r.  - 1\ .  In  this  form,  we  can  formulate  a  linear 
weighted  least-squares  problem  in  which  the  sensors  are  placed  in  a  line.  By  substituting 
r  ={rn+  q )  in  (64),  we  get 

ru\  +  lri/\  +  ri  =  Ki  -  2xixt  -  ^y,  +  Kt  ■  (65) 

By  co-locating  the  origin  with  the  first  sensor  and  making  it  the  reference  sensor,  we  get 
rt  =  /] ,  and  Kt-Kx,  then  (65)  can  be  expressed  as 

ri  +  2r</\  =  ~2xi,ixt  -  2T/,iT,  +  Kt  -  Kx .  (66) 

Finally,  with  all  the  sensors  in  a  line,  we  replace  ~2xi  Xxt -2yiXyt  by  -2x  ,  (x;  +«,>', ) 
where  at  is  some  constant.  With  this  substitution,  (66)  becomes 

r2i  - K,  +KX=  -2x.  l  (xt  +  atyt )  -  2 ri  Xrx .  (67) 
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With  (67)  now  linear  in  rx  and  [xt+aty^,  a  weighted  least-squares  solution  can  be 
obtained  as  [5],  [79] 


where 
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and  C£  is  the  inverse  covariance  matrix  of  the  TDOA  measurement  vector  from  (62). 
When  all  the  sensors  are  located  on  the  x-axis,  vn=0for  i  =  2,...,M,  resulting  in 
at  -  0.  The  solution  is  then  expressed  as 

=  (G,rc;'G,pG,rc;V  (69) 

Using  (69),  we  calculate  the  v-coordinate  estimate  using  the  expression 

<7°) 

Using  the  initial  estimate  ,  we  see  that  the  network  reorients  itself  into  a 

desired  ASLA  formation,  as  seen  in  Figure  18  (b).  In  this  configuration,  the  effect  of 
GDOP  is  minimized,  thus  minimizing  the  estimate’s  CRLB  [34].  To  take  full  advantage 
of  the  new  formation  and  its  minimized  CRLB,  we  propose  the  use  of  a  maximum- 
likelihood  estimator  introduced  in  Chapter  II.  This  estimator  is  unbiased  and 
asymptotically  efficient,  i.e.,  it  achieves  the  CRLB  [59]. 
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2.  Hyperbolic  Localization  via  Maximum-Likelihood  Estimation 

The  estimator  described  in  this  section  is  derived  from  the  generic  TDOA 
estimator  found  in  [52]  but  here  is  adapted  for  use  with  an  ASLA  formation  in  polar 
coordinates. 

Introduced  in  Chapter  II,  the  maximum-likelihood  estimator  has  been  shown  to  be 
unbiased  and  asymptotically  efficient  [59].  Its  extension  to  hyperbolic  localization  is 
done  in  a  similar  fashion  to  that  of  the  weighted  least-squares  estimator.  We  rewrite  (13) 
to  account  for  time-of-arrival  measurements  as 

t  =  t0lM+D/c  +  £  (71) 


where  D  is  an  Mxl vector  containing  the  range  from  the  emitter  to  each  node,  t0  is  the 
time  of  signal  emission,  1M  is  an  Mxl  vector  of  ones,  and  s  is  an  Mxl  vector 
containing  the  TOA  measurement  errors.  To  convert  the  TOAs  to  TDOAs,  we  subtract 
the  /lh  TOA  measurement  from  the  first  TOA  measurement  to  eliminate t0 .  This  results  in 


h  f  Di}/c  +  sl  £j . 


To  put  it  in  matrix  fonn,  we  multiply  (71)  by  the  (M-l)xAf  matrix  [52] 
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to  get 
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where  yTD  =  GMt  is  an  ( M  -  I )  x  1  vector  containing  all  the  resulting  time-difference-of- 
arrival  measurements. 

From  (18),  we  set  (or  =  6jwls  ,  and  the  maximum-likelihood  estimator’s  solution  in 
polar  coordinates  can  be  expressed  as  [52] 

=  ">wLs+c(HTGTMCt'GuHy  HTGTuC~'(yTD-(GMD„)/c),  (75) 

where  (bWLS  is  the  weighted  least-squares  estimate  in  polar  coordinates,  D0  is  the 
(M-l)xl  vector  containing  the  distances  between  the  reference  point  and  the  sensor 
nodes  and  H  is  the  (M  -  l)x  2  matrix  expressed  as 
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m  =  ,  (76) 


with  rm  and  9m  the  mth  node’s  range  and  bearing,  respectively,  r  and  Qr  are  the 
reference  point’s  range  and  bearing,  respectively,  D0m  the  mth  element  of  D0,  and  the 
product  GmH  the  hyperbolic  version  of  (17). 


B.  LOCALIZATION  PERFORMANCE 

A  series  of  Monte  Carlo  simulations  to  support  proposed  localization  technique  is 
presented  in  this  section.  The  results  of  the  simulation  are  based  on  10,000  trials  each. 
For  each  simulation,  the  ASLA  configuration  is  arranged  in  accordance  with  Figure  18, 
with  ra  =  200  m  and  N s  =  5.  Note,  since  TDOA  measurements  are  typically  on  the  order 

of  nanoseconds,  for  convenience,  we  use  range  difference-of-arrival  (RDOA),  i.e.,  the 
TDOA  measurement  weighted  by  the  signal  propagation  speed  c ,  which  is  typically  on 
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the  order  of  meters.  This  substitution  is  used  only  as  a  means  to  simplify  the  simulation 
results  and  does  not  affect  the  theory  in  any  way. 

In  Figure  19,  the  root  mean-square  error  of  location  estimate  i;loc  for  both  the 
initial  and  refined  localization  are  compared  against  their  corresponding  CRLB  as  a 
function  of  RDOA  noise  variance  <j\d  =  c1  a\  .  In  this  simulation,  Ns  =  5 ,  the  emitter  is 

located  at  (ot  =  [2500  m,  45  deg]T ,  r  =  200  m ,  and  each  data  point  is  the  result  of 

10,000  trials.  From  the  results,  we  see  that  the  4/oc  of  both  the  initial  and  refined 

estimates  increase  rapidly  with  noise.  As  expected,  the  £hc  of  the  refined  estimate 

outperforms  the  initial  estimate  at  all  values  of  noise  simulated.  Also,  we  can  see  that  the 
refined  estimator  using  a  maximum-likelihood  estimate  approaches  the  CRLB  up  to 
approximately  0  dB.  It  is  important  to  note  that  the  maximum-likelihood  estimate  is 
asymptotically  efficient  [15],  but  here,  using  the  ASLA  formation,  it  approaches  the 
CRLB  in  one  iteration. 


Figure  19.  Root  mean-square  error  %loc  of  the  initial  (black)  and  refined  (red) 
localizations  versus  RDOA  noise.  Result  based  on  10,000  trials,  Ns  =  5  , 
ra  =  200  m,  and  (ot  =  [2500  m,  45  deg]T . 
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A  plot  of  the  angular  root  mean-square  error  E,0  for  both  the  initial  and  refined 
localization  as  a  function  of  cr^,  is  provided  in  Figure  20.  In  this  simulation,  Ns  =5 ,  the 
emitter  is  located  at  cot  =  [2500  m,  45  deg]T ,  ra  =  200  m ,  and  each  data  point  is  the  result 
of  10,000  trials.  The  v-axis  is  %0  in  radians,  and  the  x-axis  is  RDOA  noise  in  dB.  From 

the  results,  we  observe  that  the  same  trends  as  in  Figure  19.  From  (56),  we  know  that  the 
ASLA  formation  yields  an  accurate  bearing  estimate.  This  is  validated  by  the  simulation 
results  with  the  refined  estimator’s  %g  = -0.25x10  ’  rad  at  -10  dB  RDOA  noise.  The 
initial  bearing  estimate  is  just  shy  of  CRLB,  while  the  refined  maximum-likelihood 
estimate  approaches  it  at  all  simulated  values.  As  with  the  <^loc  estimate,  the  maximum- 
likelihood  estimator  also  approaches  the  CRLB  in  one  iteration. 


Figure  20.  RMSE  of  the  initial  (black)  and  refined  (red)  localizations  versus 
RDOA  noise.  Result  based  on  10,000  trials,  Ns  =  5  ,  ra  —  200  m,  and 

cot  =  [2500  m,  45  deg]T . 
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We  plot  both  the  initial  and  refined  localization’s  %loc  values  in  Figure  21  as  a 
function  of  sensors  per  stack  Ns .  Each  estimator’s  CRLB  is  also  included  for 
comparison.  In  this  simulation,  cr^,  =  0  dB,  the  emitter  is  located  at 
cot  =[2000  m,  30  deg]T,  ra  =200  m,  and  each  data  point  is  the  result  of  10,000  trials. 
From  the  results,  we  see  that  the  %loc  of  the  initial  estimate  approaches  the  CRLB  at 
Ns  =2  but  quickly  diverges  from  it.  In  contrast,  the  %loc  of  the  refined  estimate  does  not 
approach  the  CRLB  until  Ns>5 .  It  is  interesting  to  note  that  an  increasing  Ns  improves 
the  performance  of  the  refined  estimate  at  a  higher  rate  than  the  initial  estimate. 


Figure  21.  Root  mean-square  error  i;loc  of  the  initial  (black)  and  refined  (red) 
localizations  versus  Ns .  Results  based  on  10,000  trials  with  ajiD  =  0  dB, 

Ns-5,  ra-  200  m ,  and cot  =  [2000  m,  30  deg]T . 

A  comparison  of  both  the  initial  and  refined  localization’s  Rvalues  is  shown  in 
Figure  22  as  a  function  of  sensors  per  stack  Ns.  In  this  simulation,  o]u:  —  0  dB.  the 
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emitter  is  located  at  cot  =  [2000  m,  30  deg]T ,  ra  -  200  m ,  and  each  data  point  is  the  result 
of  10,000  trials.  From  the  results,  we  see  that  the  of  the  initial  estimate  is  always  just 
shy  of  the  CRLB  for  all  values  of  /V  simulated.  Conversely,  the  %g  of  the  refined 
estimate  approaches  the  CRLB  for  all  values  of  Ns  simulated.  These  results  suggest  that 
the  refined  maximum-likelihood  estimator  is  a  good  candidate  for  bearing  estimates 
regardless  of  the  number  of  sensors  employed. 


Figure  22.  Root  mean-square  error  E,g  of  the  initial  (black)  and  refined  (red) 
localizations  versus  Ns .  Results  based  on  10,000  trials  wither^,  =  0  dB, 

Ns=5,  ra  =  200  m ,  and  cot  =  [2000  m,  30  deg]T . 

In  Figure  23,  both  the  initial  and  refined  localization’s  %loc  are  plotted  as  a 
function  of  the  emitter’s  bearing.  In  this  simulation,  =  0  dB,  the  emitter  is  located  at 
cot  =  [2000  m,  ,  r  =  200  m ,  and  each  data  point  is  the  result  of  10,000  trials.  From 

the  results,  it  is  clear  that  the  initial  and  refined  estimator’s  performance  is  dependent  on 
the  target’s  bearing  (relative  to  the  initial  array  configuration).  The  initial  estimate 
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follows  the  CRLB  as  seen  in  (55),  while  the  refined  estimate  improves  the  %loc  within  the 
range  of  -40  <Ot<  40  deg .  From  (55),  we  see  that  <;loc  approaches  infinity  as  6t 

approaches  ±  90  deg.  This  suggests  that  the  initial  estimate  requires  some  knowledge  of 
the  target  bearing  a  priori  to  insure  a  sufficiently  accurate  initialization  for  the  refined 
estimate.  In  other  words,  if  the  emitter  bearing  is  greater  than  ±  40  deg,  the  initial 
estimate  is  insufficient  in  setting  up  the  follow-on  refined  estimate.  If  this  is  not  possible, 
the  initial  estimate  can  be  obtained  using  estimators  that  are  independent  of  6t  such  as 
the  ones  in  [5]  and  [80].  Note  that  these  estimators  require  a  random  positioning  of  the 
sensor  nodes  and  are  not  applicable  to  an  ASLA  configuration. 
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Figure  23.  Root  mean-square  error  %loc  of  the  initial  (o)  and  refined  (*) 
localizations  versus  0t .  Results  based  on  10,000  trials  with  <y2m  -  0  dB, 
Ns=  5 ,  ra=  200  m ,  and  rt  -  2000  m . 


1,  Performance  of  the  Localization  Scheme  in  the  Presence  of  Uniform 
Position  Errors 

This  research  envisions  the  use  of  multirotor  UAVs  to  realize  an  elevated,  mobile 
WSN.  As  such,  it  is  necessary  to  address  the  position  errors  induced  by  the  UAV’s 
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unsteady  station-keeping  operations.  This  unsteadiness  is  mainly  due  to  wind  [46]  and 
causes  nodes  to  momentarily  deviate  from  their  given  positions.  In  this  section,  we 
model  these  position  errors  as  two  unifonnly  distributed  random  variables  Sx  and  8  , 

each  bounded  by  8  .  The  true  distribution  of  position  errors  for  any  given  UAV  is 

determined  by  many  factors  such  as  GPS  accuracy,  wind  effect,  and  flight  controller  [38]. 
In  the  absence  of  such  knowledge,  we  resort  to  a  basic  distribution  in  order  to  complete 
our  analysis.  The  effects  of  these  position  errors  are  illustrated  in  Figure  24  through 
Figure  27. 

A  comparison  of  the  initial  localization’s  root  mean  squared  location  error  £loc 
and  bearing  error  E,e  with  and  without  position  errors  is  shown  in  Figure  24  and 
Figure  25  as  a  function  of  the  maximum  position  error  8p  (see  (53)).  In  this  simulation, 

a\D  =  0  dB,  the  emitter  is  located  at  a>t  =  [2000  m,  30  deg]T ,  ra  -  200  m ,  and  each  data 
point  is  the  result  of  10,000  trials.  From  the  plots,  we  see  that  the  root  mean  squared 
location  error  monotonically  increases  with  8p ,  both  in  location  and  bearing  estimates. 


Figure  24.  Initial  localization’s  root  mean-square  error  i;loc  with  no  position  errors 
(black)  and  with  position  errors  (red)  versus  maximum  uniform  position 
error  Sp .  Results  based  on  10,000  trials  with  crj^  =  0  dB,  Ns  =  5 , 

ra  =  200  m ,  and  cot  =  [2000  m,  30  deg]T . 
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Figure  25.  Initial  localization’s  root  mean-square  error  ^  with  no  position  errors 
(black)  and  with  position  errors  (red)  versus  maximum  uniform  position 
error  5p .  Results  based  on  10,000  trials  with  cr2^  =  0  dB,  Ns=  5, 

r  =  200  m ,  and  cot  =  [2000  m,  30  deg]T . 

Similar  to  Figure  24  and  Figure  25,  the  performance  of  the  refined  localization  as 
a  function  of  8  is  shown  in  Figure  26  and  Figure  27.  In  this  simulation,  crRj)  —  0  dB,  the 

emitter  is  located  at  cot  =  [2000  m,  30  deg]  ,  ra  =  200  m ,  and  each  data  point  is  the  result 
of  10,000  trials.  From  the  results,  we  see  that  the  refined  maximum-likelihood  estimate 
is  more  resilient  against  position  error  with  an  almost  linear  response  to  Sp  .  Specifically, 

we  see  in  Figure  26  that  for  8  —  1  m,  <^loc  only  increases  by  approximately  50  m. 

Overall,  these  results  suggest  that  position  errors  play  a  significant  role  in  the  accuracy  of 
localization  using  a  mobile  network.  In  the  next  section,  we  address  these  errors  using  a 
measurement  outlier  rejection  process. 
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Figure  26.  Refined  localization’s  root  mean-square  error  £loc  with  no  position 
errors  (black)  and  with  position  errors  (red)  versus  maximum  uniform 
position  error  8  .  Results  based  on  10,000  trials  with  a ^  =  0  dB,  Ns=5, 

ra  =  200  m ,  and  cot  =  [2000  m,  30  deg]T . 


Sp(m) 


Figure  27.  Refined  localization’s  root  mean-square  error  %6  with  no  position  errors 
(black)  and  with  position  errors  (red)  versus  maximum  uniform  position 
error  Sp .  Results  based  on  10,000  trials  with  =  0  dB,  Ns  =  5 , 

ra  =  200  m ,  and  cot  =  [2000  m,  30  deg]T . 
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C.  ROBUST  SIGNAL  LOCALIZATION  VIA  MEASUREMENT  OUTLIER 

REJECTION 

With  the  goal  of  improving  the  robustness  of  the  source  localization  technique, 
we  propose  a  measurement  outlier  rejection  process  to  improve  localization  noise 
performance.  This  includes  resilience  to  increased  signal  noise,  sensor  position  errors, 
and  measurement  outliers.  These  improvements  leverage  the  existing  outlier  detection 
techniques  to  identify  and  remove  specious  measurement  data,  i.e.,  measurements  that 
seem  benign  but  in  actuality  degrade  the  estimate’s  accuracy.  The  basic  premise  is  to 
determine  the  effect  each  measurement  has  on  the  end  solution  as  a  means  to  identify  and 
remove  measurement  outliers.  In  this  research,  we  define  a  measurement  outlier  as  a 
measurement  taken  with  a  variance  greater  than  a  specified  multiple  of  the  expected 
variance. 

1.  Robust  Hyperbolic  Localization  via  Measurement  Outlier  Rejection 

When  trying  to  detect  the  presence  of  outliers  in  a  data  set  containing  location 
estimates,  the  simplest  methodology  is  to  calculate  each  sample’s  Euclidean  distance 
from  the  sample  mean.  The  problem  with  this  approach  is  that  outliers  heavily  influence 
the  sample  mean,  making  it  a  poor  discriminator  of  outliers  [18].  Another  approach  is  to 
use  a  statistical  hypothesis  test  to  determine  the  likelihood  that  the  data  set  is  derived 
from  a  single  distribution  [55];  however,  this  approach  is  more  suited  to  problems  with 
larger  data  sets  (>  30  samples  per  stack)  [55],  which  is  not  appropriate  for  our  scheme. 

a.  Measurement  Outlier  Rejection  via  the  Mahalanobis  Distance 

Since  the  localization  technique  is  focused  on  a  small  sample  size,  the  use  of 
hypothesis  testing  to  identify  outliers  is  limited.  An  alternative  approach  is  the 
Mahalanobis  distance.  This  approach  performs  well  with  small  sample  sizes  and  is  less 
affected  by  the  presence  of  outliers  [18].  In  contrast  to  the  Euclidean  distance,  the 
Mahalanobis  distance  calculates  the  distance  from  the  mean  but  does  so  in  a  manner  that 
accounts  for  the  samples  distribution,  i.e.,  its  covariance  matrix.  Given  Nsc  location 
estimates,  the  corresponding  Mahalanobis  distance  is  defined  as  [18] 
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dM,i  =  (^i-E{va})Ca-1  (v,  -E{Vco})T , 


(77) 


where  the  vwis  the  Nsc  x 2  matrix  containing  all  the  location  estimates,  v.  is  the  1x2 

vector  containing  the  7th  location  estimate,  and  CJs  the  2x2  covariance  matrix  of  vm. 

In  the  proposed  rejection  process,  we  use  the  Mahalanobis  distance  of  the  weighted  least- 
squares  solution  to  determine  the  best  and  worst  candidate  measurements  for  use  in  the 
final  location  estimate.  In  this  approach,  we  use  single  case  diagnostics  [18]  to  derive 
multiple  weighted  least-squares  estimates  from  one  set  of  TDOA  measurements.  For 
ASLA  configured  network,  the  number  of  single  cases  for  a  given  signal  measurement  is 

NX=M- 1. 

A  classic  tool  for  the  detection  of  statistical  outliers  in  least-square  estimation 
problems,  single  case  diagnostics  provide  a  methodology  to  single  out  measurements  that 
lead  to  strong  deviations  in  the  end  solution  [18].  In  this  approach,  a  single  measurement 
from  a  set  is  removed,  and  the  corresponding  weighted  least-squares  solution  is 
calculated.  This  process  is  repeated  until  each  measurement  is  treated  in  the  same 
fashion.  The  Mahalanobis  distance  of  each  solution  is  then  calculated  using  the  mean 
and  covariance  of  the  entire  Nsc  set  of  solutions.  The  solutions  with  high  Mahalanobis 

distance  values  suggest  that  an  influential  measurement  was  removed.  Such 
measurements  are  termed  leverage  points  [18].  Conversely,  if  the  Mahalanobis  distance 
value  was  low,  then  the  measurement  in  question  is  considered  valid. 

b.  Distribution  of  the  Mahalanobis  Distance 

With  a  measure  of  solution  divergence  now  obtained,  the  next  step  is  to  be  able  to 
discern  which  values  of  this  divergence  represent  outliers  and  which  represent  valid 
measurements.  This  process  amounts  to  setting  a  Mahalanobis  distance  value  outlier 
threshold.  To  calculate  a  suitable  threshold  value,  we  consider  the  distribution  of  the 
Mahalanobis  distance  values  with  and  without  outliers.  Given  a  set  of  location  estimates 
that  are  Gaussian  with  a  known  mean  and  covariance  matrix,  the  distribution  of  the 
Mahalanobis  distance  values  follow  a  Chi-squared  distribution  [81].  The  Chi-squared 
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distribution  is  the  resulting  distribution  for  the  sum  of  multiple  squared  Gaussian  random 
variables.  The  Chi-squared  distribution  can  be  expressed  by  its  degrees  of  freedom  p, 
which  is  the  number  of  squared  Gaussian  random  variables  present  in  the  sum.  The  Chi- 
squared  distribution  of  the  Mahalanobis  distance  is  expressed  as  [55] 

fnM  K/)=  2pI2T(PI2)  6XP ' l2)idM  )/,/2~‘  for  0  <  dM  <  oo ,  (78) 


where  T  is  the  gamma  function  defined  as 


r(r)  =  j xT  le'dx  for  r  >  0 . 

o 

The  cumulative  distribution  function  is  defined  as 

df  x^-2)/2 e~x'2 

and  the  corresponding  means  and  variance  are  expressed  as  [55] 


E{dii}  =  P> 

a2d  =2  p. 

aM 


(79) 


(80) 


(81) 


Given  that  our  scenario  considers  a  small  number  of  outliers,  we  hypothesize  that  the 
distribution  of  the  Mahalanobis  distance  with  and  without  outliers  is  significantly  similar. 

To  validate  this  hypothesis,  a  comparison  of  Mahalanobis  distance  values  and 
their  estimated  distribution  are  shown  in  Figure  28.  In  this  simulation,  using  (69),  the 
histogram  of  100  Mahalanobis  distance  values  with  and  without  outliers  is  shown.  The 
Chi-squared  probability  density  functions  were  plotted  using  the  sample  mean.  In  the 
outlier  case,  there  are  15  outliers  present  with  a  variance  five  times  larger  than  the 
expected  variance.  In  both  cases,  the  estimated  degrees  of  freedom  was  yz  =  1.98.  These 
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results  help  support  the  claim  that  the  estimated  probability  density  function  of  the 
Mahalanobis  distance  approximates  the  Chi-squared  with  or  without  a  small  amount  of 
outliers  present. 


Mahanobis  distance 

a)  Histogram  and  PDF  of  Mahalanobis  distance  values  without  outliers 
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b)  Histogram  and  PDF  of  Mahalanobis  distance  values  with  outliers 
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Figure  28.  Comparison  of  histogram  and  Chi-squared  probability  density  function 

of  Mahalanobis  distances. 


c.  Outlier  Rejection  Threshold 

Using  the  previous  results  to  support  our  hypothesis,  we  next  determine  the 
outlier  rejection  threshold.  A  standard  criterion  for  outlier  rejection  thresholds  is  the 
Tukey  criteria  for  anomalies  [82].  This  criterion  states  that  in  a  given  range  of  data 
points,  values  at  the  extreme  ends  of  the  data  set  are  considered  outliers.  To  quantify  the 
threshold,  Tukey  [82]  proposed  the  following  criterion.  Taking  the  difference  between 
the  25th  and  75th  percentile  of  a  data  set  as  the  interquartile  range  [82] 


Q  m  Q?5%  Q  25% 
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(82) 


we  consider  any  data  point  greater  than  Q75o/o+  (3/2)  Q//(  as  an  outlier  [82],  It  is  important 

to  note  that  this  criterion  makes  no  consideration  on  the  data’s  distribution.  To 
incorporate  knowledge  of  the  data  set  distribution,  we  propose  the  use  of  Horn’s 
algorithm  [83], 

Horn’s  algorithm  begins  with  the  estimation  of  a  data  set’s  distribution.  Next,  the 
Tukey  criterion  threshold  is  detennined  using  the  estimated  cumulative  distribution 
function.  Then,  all  data  points  greater  than  the  calculated  threshold  are  removed  from  the 
data  set  [83].  Overall,  assuming  the  distribution  of  the  Mahalanobis  distance  values  is 
Chi-squared,  the  interquartile  range  is  detennined  by  [83] 

Qffl  =  F1  (0.75 1  p)  -  F1  (0.25  |p),  (83) 

where  the  F  1  is  the  inverse  Chi-squared  cumulative  distribution  function.  The  outlier 
threshold  using  Horn’s  algorithm  can  then  be  expressed  as  [82]  [83] 

T„m  =  F-'(0.75|p)  +  (3/2)Qa.  (84) 

Using  the  outlier  rejection  process  in  conjunction  with  our  proposed  localization, 
we  wish  to  minimize  the  effects  of  TDOA  measurement  noise,  position  errors,  and 
measurement  outliers.  We  term  an  estimator  using  this  process  a  robust  estimator. 

2.  Performance  of  the  Robust  Localization  in  the  Presence  of 
Measurement  Outliers 

To  support  the  proposed  outlier  rejection  process,  we  examine  the  results  of  four 
Monte  Carlo  simulations.  Each  simulation  evaluates  the  robust  localization  estimator’s 
performance  in  the  presence  of  random  measurement  outliers.  Each  simulation  uses  a  15 
node  ASLA  configured  network  in  which  the  emitter  is  located  at 
cot  =[2000  m,  30  deg]T. 
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A  comparison  of  the  robust  and  regular  weighted  least-squares  estimator’s  root 
mean-square  error  %loc  as  a  function  of  <y2t/)  =  c2 a\  is  shown  in  Figure  29.  In  this 
simulation,  there  are  no  measurement  outliers  present,  N s  =  5 ,  the  emitter  is  located  at 
cot  =[2000  m,  30  deg]T,  ra  -  200  m,  and  each  data  point  is  the  result  of  10,000  trials. 
The  results  show  that  the  robust  estimator’s  £hcis  only  slightly  higher  than  the  regular 

estimator  at  low  noise  levels.  In  contrast,  the  robust  estimator  outperforms  the  regular 
estimator  at  high  noise  levels.  Overall,  this  result  validates  the  effectiveness  of  the 
rejection  process  against  high  measurement  noise. 


Figure  29.  Root  mean-square  error  %loc  of  the  weighted  least-squares  (red)  and 
robust  weighted  least-squares  (blue)  versus  RDOA  noise  variance.  No 
outliers.  Results  based  on  10,000  trials  with  Ns  =  5 ,  ra=  200  m ,  and 

cot  =  [2000  m,  30  deg]T . 

In  Figure  30,  the  robust  weighted  least-squares  estimator’s  %loc  is  compared 

against  its  corresponding  CRLB  at  increasing  RDOA  noise  values.  In  this  simulation, 
there  are  two  random  measurement  outliers  with  a  RDOA  noise  variance  that  is  five 
times  greater  than  the  remaining  set  of  measurements.  Also,  Ns  -  5 ,  the  emitter  is 
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located  at  cot  =  [2000  m,  30  deg]T,  ra  -  200  m ,  and  each  data  point  is  the  result  of  10,000 

trials.  The  v-axis  is  %loc  in  meters,  and  the  x-axis  is  RDOA  noise  variance  in  dB.  The 

results  show  that  the  effects  of  the  outliers  are  almost  negligible  at  the  lower  levels  of 
RDOA  noise  (<  -5  dB).  The  robust  weighted  least-squares  estimator  outperforms  the 
weighted  least-squares  estimate  at  all  simulated  noise  levels,  with  significant 
improvements  occurring  at  RDOA  noise  level  above  -2  dB  .  Of  interest  is  the  robust 
estimator’s  performance  above  6  dB,  where  it  again  outperforms  the  no  outlier  case.  This 
result  validates  the  effectiveness  of  the  proposed  outlier  rejection  process  for  a  wide 
range  of  measurement  noise,  with  or  without  outliers  present. 


Figure  30.  Root  mean-square  error  i;loc  of  the  weighted  least-squares  estimator  with 
outliers  (red),  weighted  least-squares  estimator  without  outliers  (black),  and 
robust  weighted  least-squares  estimator  with  outliers  (blue)  versus  RDOA 
noise  variance.  Results  based  on  10,000  trials  with  Ns  =  5 ,  ra=  200  m , 

and  a>t  =  [2000  m,  30  deg]T .  In  the  outlier  cases  there  were  two  random 
outliers  present  at  a  RDOA  noise  variance  five  times  <j\D . 
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The  results  of  two  Monte  Carlo  simulations  illustrating  the  robust  weighted  least- 
squares  estimator’s  performance  as  a  function  of  the  number  of  random  outliers  present 
and  their  variance  are  shown  in  Figure  3 1  and  Figure  32. 

A  plot  of  the  robust  and  regular  weighted  least-squares  estimator  £/oc  is  shown  in 
Figure  31  as  a  function  of  the  number  of  random  outliers  present.  In  this  simulation, 
Ns  =  5 ,  the  emitter  is  located  at  cot  =  [2000  m,  30  deg]T ,  ra  =  200  m ,  c>)in  =  0  dB ,  the 
outlier  variance  is  five  times  larger  than  the  expected  variance,  and  each  data  point  is  the 
result  of  10,000  trials.  From  the  results,  we  see  that  the  increase  in  E,loc  of  both  the 

regular  and  robust  estimate  is  almost  linear,  with  the  robust  estimator  having  a  slower 
rate  of  increase.  This  result  suggests  that  the  performance  increase  of  the  robust 
estimator  is  a  function  of  the  number  of  outliers  present.  As  the  number  of  outliers 
increases,  the  more  of  them  are  rejected. 


Figure  31.  Root  mean-square  error  %loc  of  the  weighted  least-squares  estimator  with 
outliers  (red),  weighted  least-squares  estimator  without  outliers  (black),  and 
robust  weighted  least-squares  estimator  with  outliers  (blue)  versus  number 
of  outliers.  Results  based  on  10,000  trials  with  Ns  -  5  ,  ra=  200  m ,  and 

cot  =  [2000  m,  30  deg]T .  In  the  outlier  cases  their  RDOA  noise  variance  is 

five  times  crL . 


65 


A  comparison  of  the  robust  and  regular  weighted  least-squares  estimator’s  %loc 

values  as  a  function  of  outlier  variance  gain  is  shown  in  Figure  32.  In  this  simulation,  the 
parameters  are  similar  to  the  ones  used  to  create  Figure  31.  From  the  results,  we  see  that 
the  increase  in  %loc  approaches  a  maximum  at  variance  multiplier  greater  than  ten.  Since 
each  outlier  is  derived  with  a  larger  variance,  they  are  easier  to  identify  and  be  removed, 
thus  we  see  a  plateauing  effect  on  the  %loc.  This  also  suggests  that  above  a  certain 

variance,  value  the  %loc  of  the  robust  estimator  has  a  stronger  dependence  on  the  number 
of  outliers  rather  than  on  their  variance. 


Oultier  Measurements  RDOA  noise  variance  (a^D) 

Figure  32.  Root  mean-square  error  %Ioc  of  the  weighted  least-squares  estimator  with 
outliers  (red),  weighted  least-squares  estimator  without  outliers  (black),  and 
robust  weighted  least-squares  estimator  with  outliers  (blue)  versus  RDOA 
noise  gain.  Results  based  on  10,000  trials  with  Ns  —  5  ,  ra  -  200  m ,  and 

a>t  =  [2000  m,  30  deg]T .  In  the  outlier  cases  there  were  2  random  outliers 

present. 

The  performance  of  the  robust  weighted  least-squares  estimator  compared  to  the 
regular  weighted  least-squares  estimator  as  a  function  of  the  emitter’s  bearing  is  shown  in 
Figure  33.  In  this  simulation,  there  are  two  outlier  measurements  present  with  a  variance 
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equal  to  five  times  a 1RD,  the  emitter  is  located  at  a>t  -  [2000  m,  where  0t  ranges 
from  -60  <  6t  <  60 ,  ra  -  200  m ,  and  each  data  point  is  the  result  of  10,000  trials.  From 
the  results,  we  see  the  robust  weighted  least-squares  estimator  performs  better  than  the 
regular  estimator  and  at  a  wider  range  of  0t .  This  suggests  that  the  robust  estimator  has  a 
smaller  dependence  on  the  emitter’s  bearing  compared  to  the  regular  weighted  least- 
squares  estimate.  This  also  validates  the  effectiveness  of  the  rejection  processes  in 
increasing  the  initial  estimate’s  accuracy  in  the  presence  of  measurement  outliers.  From 
Figure  23,  we  see  that  the  refined  estimator’s  performance  is  dependent  on  the  initial 
estimator’s  perfonnance;  therefore,  any  improvement  to  the  initial  estimate’s  accuracy 
benefits  the  overall  localization  perfonnance. 


300 

280 

260 

240 

220 

—  200 
8  180 
160 
140 
120 
100 
80 


Figure  33.  Root  mean-square  error  <;hc  of  the  weighted  least-squares  estimator 
(circle)  and  robust  weighted  least-squares  estimator  (diamond)  versus  0t . 
Results  based  on  3,000  trials  with  cr^,  =  0  dB,  Ns  —  5  ,  r  =  200  m ,  and 
rt  =  2000  m .  In  both  cases  there  were  2  random  outliers  present  at  a  RDOA 

noise  variance  five  times  aln . 
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D.  PERFORMANCE  OF  THE  ROBUST  LOCALIZATION  IN  THE 
PRESENCES  OF  MEASUREMENT  OUTLIERS  AND  UNIFORM 
POSITION  ERRORS 

In  this  section,  we  examine  at  the  performance  of  the  proposed  robust  weighted 
least-squares  localization  in  the  presence  of  measurement  outliers  and  uniform  position 
errors.  To  assess  the  performance  of  the  localization  with  position  errors  we  evaluated 
two  Monte  Carlo  simulations,  with  and  without  measurement  outliers  present.  Each 
simulation  uses  a  15  node  ASLA  configured  network  in  which  the  emitter  is  located  at 
a>t  =  [2000  m,  30  deg]T ,  r  =  200  m ,  0)U)  =  0  dB, ,  and  all  data  points  are  the  result  of 
10,000  trials. 


A  comparison  of  the  robust  and  regular  weighted  least-squares  estimator’s  %loc  as 
a  function  of  maximum  uniform  position  error  8  is  shown  in  Figure  34.  In  this  case, 
there  are  no  measurement  outliers  present.  From  the  results,  we  see  that  the  robust 
estimator  does  not  improve  performance  over  the  regular  estimator  until  a 8  value  of  0.4 

meters.  Overall,  this  would  suggest  that  the  outlier  rejection  process  is  not  beneficial  in 
the  small  position  error  case  given  that  there  are  no  outliers  present. 


Figure  34.  Root  mean-square  error  %loc  of  the  weighted  least-squares  estimate  (red) 
and  robust  weighted  least-squares  estimate  (blue)  versus  8  .  For 
comparison  the  no  error  case  is  shown  in  black.  Results  based  on  10,000 
trials  wither^  =  0  dB,  Ns  =  5 ,  ra  —  200  m ,  and  (ot  =  [2000  m,  30  deg]T . 
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Now,  with  measurement  outliers  present,  a  comparison  of  the  robust  and  regular 
weighted  least-squares  estimator’s  %loc  is  shown  in  Figure  35.  In  contrast  to  Figure  34, 

the  robust  estimator  shows  a  semi-fixed  amount  improvement  over  the  regular  estimator 
at  all  values  of  position  errors  simulated.  This  result  suggests  that  the  robust  estimator  is 
effective  against  measurement  outliers  but  only  marginally  effective  in  the  presence  of 
positon  errors. 


6p(m) 

Figure  35.  Root  mean-square  error  %loc  of  the  weighted  least-squares  estimate  (red) 
and  robust  weighted  least-squares  estimate  (blue)  versus  8  .  For 
comparison  the  no  error  case  is  shown  in  black.  Results  based  on  10,000 
trials  wither^,  =  0  dB,  Ns  =  5 ,  ra  =  200  m ,  and  cot  =  [2000  m,  30  deg]T . 

In  this  chapter,  we  developed  a  localization  technique  that  is  capable  of 
approaching  the  CRLB.  This  technique  consists  of  two  sequential  location  estimators,  a 
weighted  least-squares  estimator  followed  by  a  maximum-likelihood  estimator  to  refine 
the  initial  estimate.  We  demonstrated  that  this  technique  is  moderately  susceptible  to 
sensor  position  errors.  Next,  we  developed  a  measurement  outlier  rejection  process  to 
compensate  for  these  position  errors  as  well  as  measurement  outliers.  The  simulation 
results  showed  that  the  outlier  rejection  technique  is  effective  against  measurement 
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outliers  and  large  position  errors.  In  the  next  chapter,  we  provide  a  detailed  description 
of  the  proposed  collaborative  signal  collection  technique. 
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V.  ROBUST  SIGNAL  COLLECTION  FROM  AN  AIRBORNE 
SEMI-STATIONARY  NETWORK 


In  this  chapter,  we  develop  the  proposed  collaborative  signal  collection  scheme. 
This  begins  with  an  analysis  of  uniform  and  Gaussian  sensor  position  errors  and  their 
effects  on  the  ASLA  formation’s  array  factor  and  main  beam  gain.  Through  this  analysis, 
we  are  able  to  derive  a  novel  signal  estimator  that  uses  sampled  data  and  knowledge  of 
the  position  error  statistics  to  compensate  for  the  sensor  position  errors.  We  further 
improve  the  collection’s  robustness  to  interfering  signals  by  adapting  virtual  filling  and 
array  tapering  techniques  to  the  proposed  scheme.  Simulations  are  used  to  support  the 
developed  theory  and  analyze  its  performance. 

A.  EFFECTS  OF  UNIFORM  POSITION  ERROR  ON  THE  ARRAY  FACTOR 

In  Chapter  IV,  we  showed  that  uniform  position  errors  negatively  affect  the 
localization  performance.  In  this  section,  we  discuss  the  effects  of  these  errors  on 
collaborative  beamforming. 

From  [63],  we  rewrite  (28)  to  represent  the  array  factor  of  a  planar  array  as 

=  (85) 

where  0t  is  the  signal  bearing ,6sa  is  the  array  steering  angle,  i.e.,  the  direction  of  the 
main  beam  gain,  at  is  the  7th  node’s  signal  phase  expressed  as 

ai  =  [5  (jqsin#,  +  y;cos6*, ) ,  (86) 

wi  is  the  7th  node’s  complex  weight  expressed  as 

=  exP  (~JP  (*,-sin0m  +  J,cos6(a )) ,  (87) 
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M  is  the  number  of  sensor  nodes,  and  /?  =  2 n  /  X  .  Incorporating  unifonn  position  errors 
Sx  and  8y  (see  (53))  into  (86)  yields 


«  =  P[(xt  +  5xi )  sin  <9,  +  (y,  +  Syi )  cos  6, ) .  (88) 

Separating  the  sensor  position  from  its  position  errors,  we  get 

As  {°SaAA>Sy)  =  ZmW/CXP  (M  )eXp(^-).  (89) 

where 

aPj  =  +  SyjCosO, )  (90) 

is  the  Ith  node’s  phase  perturbation  due  to  position  errors.  To  steer  the  array  to  the 
intended  signal’s  direction  0t,  we  set  6sa  equal  to  6 t .  As  a  result,  vycxp (  jai )  =  1 ,  and 
the  array  factor  becomes  the  expression  for  the  main  beam  gain 


y|t 

0sa  =0,  =  constant  1  ~ 1 


M  jf3(8x  is\nOt  +dyicosOt ) 


(91) 


where  we  assume  that  0t  is  known  and  treat  it  as  a  constant.  By  letting 


jplS.iSind.+S v.cosft)  ,,  i  i 

>Jei=e  ,  we  can  express  the  mean  ot  the  phase  error  per  node  as 


{Sr)dS^Srl . 


(92) 


By  substituting /A  (Sx)  =  fA  [Sy )  =  1  / [25 p  j  into  (92),  we  obtain 
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magnitude  \/i 


nsi 


sin  (^SpPsmOt  )  sin  p/3  cos  &t  j 
(Spfif  sinO'COsO, 


and 


angle  \u  =0, 


(93) 


(94) 


The  variance  of  //, ;  is  then  obtained  from 


^  = 


(95) 


where  the  second  moment  is  j/^.j  is  expressed  as 


£>„  <•<?„ 


,-sin^  +^,  j-cos^ ) 


f  .  \ 


5p  C^p  j /?( 8X , sin 0t  +Sy  /Cos Gt  )  -jfi(dx  is\nOt+dy  ,cos^) 


v2^y 


f  ,  A 


v2A 


Cl 


dS  dS  , 

x,i  y,i  ’ 


(96) 


=  1. 


Using  (96),  we  rewrite  (95)  as 


(97) 


Assuming  that  8X  .  and  8  .  are  independent  and  identically  distributed  random  variables, 
we  have  that  rjEi  for  i  =  are  also  independent  and  identically  distributed.  Given 
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that  GMB[e„dx,dy)  is  the  sum  of  these  M  random  variables,  we  express  its  mean  and 
variance  as 

M 

»cm  =IX,  =M^ni:  (98) 

/= 1 

and 

M 

(*>> 

i= 1 

1,  Effects  of  Uniform  Position  Error  in  the  Special  Case  when  the  ASLA 
Formation  Is  Normal  to  the  Target  Source 

In  the  special  case  where  the  ASLA  array  has  been  reoriented  (see  Figure  18  (b)), 

the  target  emitter  is  normal  to  the  array  making  6t  -  6sa  —  0  deg.  In  this  scenario,  the 

expressions  for  the  mean  of  the  main  beam  gain  is  simplified  such  that  (91)  is  now  only  a 

function  of  5  and  is  written  as 

y 


M  JPSy. 


<S»K)  =  L",e 


(100) 


By  letting  fj '  ;  =  e  'pSyi ,  we  can  express  the  mean  of  the  phase  error  per  node  as 


Mn 


(101) 


By  substituting  /A  (S, )  =  1  /  {28 p  j  into  (101),  we  obtain 


magnitude  \fi. 


sin 


8J 


(102) 


angle  =0. 


(103) 


Along  the  lines  of  (98)  and  (99),  we  express  mean  and  variance  of  the  main  beam  gain  as 


74 


(104) 


M 

A.  =2X  =mK 

1=1 

and 

M 

<105> 

1=1 

2.  Simulation  Results  on  the  Effects  of  Uniform  Position  Errors  on  the 
Array  Factor 

To  support  these  findings,  three  simulations  were  conducted  using  an  N s  =  5 
ASLA  configured  network  with  a  total  of  M  - 15  nodes  each  with  unifonn  position 
errors.  For  these  simulations,  A  =  1 ,  the  source  emitter  is  located  at  a  bearing  of 
6t  =  45  deg ,  ra  =  200  m ,  and  all  data  points  are  the  result  of  1,000  trials. 

A  comparison  between  the  simulated  and  theoretical  values  of  jum  and  as  a 

function  of  Sp  are  shown  in  Figure  36  and  Figure  37.  From  the  results,  we  see  that  the 

simulated  values  closely  adhere  to  the  theoretical  values,  thus  validating  the  derived 
expressions. 


Figure  36.  Mean  normalized  main  beam  gain  of  an  Ns  =  5  ASLA  versus  8  . 
Results  based  on  1,000  trials  withr  =  200  m ,  and  0t  =  45  deg  . 
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Figure  37.  Variance  of  the  normalized  main  beam  gain  for  Ns  =  5  ASLA  versus 
Sp  .  Results  based  on  1,000  trials  withra  =  200  m ,  and  9t  =  45  deg  . 

We  now  evaluate  the  effects  of  unifonn  position  errors  on  the  array  factor.  The 
average  normalized  array  factor  error  (//  .  (refer  to  (61))  as  a  function  of  maximum 

uniform  position  error Sp  is  shown  in  Figure  38.  In  this  simulation,  Ns=5,  each  data 
point  is  the  result  of  1,000  trials,  2  =  1,  ra-  200  m ,  and  6t  =  45  deg  .  From  the  results, 
we  can  see  that  the  average  y/E  rapidly  increases  with  8  ,  which  means  that  even  minor 

position  (Sp  <1  mj  errors  can  have  a  significant  effect  on  the  array  factor.  This  suggests 

that  position  errors  of  this  nature  if  not  compensated  for  will  negate  any  benefits  that 
collaborative  beamforming  can  provide.  This  result  is  a  key  motivator  for  this  research, 
as  it  is  the  objective  of  this  dissertation  to  deliver  an  effective  signal  collection  scheme 
(see  Chapter  V  section  C)  in  the  presence  of  such  errors. 
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Figure  38.  Average  y/£  for  an  Ns  =  5  ASLA  versus  8  .  Results  based  on  1,000 
trials  with  X  =  1,  ra  =  200  m ,  and  6t  =  45  deg  . 

B.  EFFECTS  OF  GAUSSIAN  POSITION  ERROR  ON  THE  ARRAY 

FACTOR 

In  the  previous  section,  we  evaluated  the  effects  of  unifonn  position  error  on  the 
array  factor  and  the  main  beam  gain  of  an  ASLA  configuration.  In  this  section,  we 
perform  a  similar  analysis  with  Gaussian  position  errors. 

Similar  to  the  unifonn  position  error  case,  we  wish  to  obtain  an  expression  for  the 
mean  and  variance  of  the  main  beam  gain;  however,  using  the  same  methodology 
becomes  intractable  for  the  Gaussian  random  variables.  As  an  alternative,  we  instead 
evaluate  the  main  beam  gain  after  the  ASLA  formation  has  been  reoriented  as  shown  in 
Figure  18  (b).  In  this  configuration,  we  assume  that  6t  -  9sa  =  0  deg.  This  assumption 
allows  us  to  study  the  effects  of  Gaussian  position  errors  independent  of  errors  in  the 
bearing  estimate  9t .  Asa  result  of  this  assumption,  the  main  beam  gain  from  (91)  can  be 
expressed  as 
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(106) 


M  jps 


G^SJ^e 


jj3Sy 

where  e  is  known  as  the  wrapped  Gaussian  random  variable  [84].  By  letting 

_  jPSyi 

,  we  can  interpret  the  main  beam  gain  as  the  sum  of  M  wrapped  Gaussian 
random  variables'4^,  each  with  a  mean  of  zero  and  variance  a]  =  ( /3a p  j  [84],  The 
probability  density  function  of  the  wrapped  Gaussian  random  variable  is  given  by  [84] 


a..K-«J)=t^7ZL.«p 


2  71(7  „ 


(wg  +  2nk^ 

2  a! 


2  A 


(107) 


The  mean  values  of  the  magnitude  and  phase,  respectively,  are  given  by  [84] 


magnitude  { juw  J  =  e  n  (108) 

angle  {//w}  =  0  .  (109) 

The  variance  is  given  as 

ol=\-e^n.  (110) 

Since  the  main  beam  gain  is  the  sum  of  M  terms,  similar  to  (98)  and  (99)  we  express  the 
mean  value  of  its  magnitude  and  phase  as 


magnitude  { juMB  }  =  Me  *  n , 

(111) 

angle}  fiMB }  =  0 

(112) 

with  a  corresponding  variance  of 

a-2  =m{  l-e'(/?t7p)2/2l 

u  MB  1V1  1  c 

V  J 

(113) 
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1.  Simulation  Results  of  the  Effects  of  Gaussian  Position  Errors  on  the 
Array  Factor 

To  support  the  theory  presented  in  the  previous  section,  three  simulations  were 
conducted  using  an  Ns  =30  ASLA  configured  network,  with  Gaussian  position  errors. 

In  these  simulations,  the  source  emitter  is  located  at  a  bearing  of  0t  =0deg,  2  =  1, 
ra  =  200  m ,  and  all  data  points  are  the  result  of  1,000  trials. 

A  plot  of  the  simulation  and  theoretical  normalized  values  of  magnitude  { /uMB  }  as 

a  function  of  cr  =(/?cr/;)  is  shown  in  Figure  39.  From  the  plot,  we  see  that  the 

simulation  results  closely  follow  the  theoretical.  The  results  indicate  that  the  main  beam 
gain  monotonically  decreases  as  cr  increases.  Similar  to  the  uniform  position  error  case, 

we  see  that  the  errors  significantly  reduce  the  effectiveness  of  collaborative 
beamforming. 


Figure  39.  Normalized  values  of  magnitude  [juMB  }  versus  position  error  variance 

cr  =  ( j3(Tp )  .  In  this  simulation  Ns  -  30 ,  6t  =  0  deg  ,2  =  1,  r  =  200  m , 
and  all  data  points  are  the  result  of  1,000  trials. 
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Similar  to  the  previous  simulation,  a  comparison  between  the  simulation  and 
theoretical  values  of  angle as  a  function  of  cr  =  )  is  shown  in  Figure  40. 

From  the  results,  we  see  that  the  simulation  results  closely  follow  the  theoretical  value, 
with  the  angle {juMB}  values  approximately  equal  to  zero  at  all  values  of  cr  simulated. 
This  suggests  that  the  effects  of  Gaussian  position  errors  have  no  effect  on  the  mean 
value  of  the  main  beam  phase  angle.  Since  position  errors  do  not  affect  the  angle  {/uMB)  , 
it  will  have  a  prominent  role  in  the  signal  estimation  described  in  the  following  section. 


(P<Jpr  (rad  ) 


Figure  40.  Values  of  angle  versus  position  error  variance  g  yr  p>  .  In  this 

simulation  Ns  =  30 ,  0t  -  0 deg ,  X  =  1,  ra=  200  m ,  and  all  data  points  are 

the  result  of  1,000  trials. 


C.  SIGNAL  COLLECTION  FROM  AIRBORNE  SYMMETRIC  LINE  ARRAY 
NETWORK  IN  THE  PRESENCE  OF  UNIFORM  POSITION  ERRORS 

In  the  previous  section,  we  showed  that  each  node’s  phase  perturbations  due  to 

position  errors  have  a  mean  value  of  zero.  Given  that  the  ASLA  configuration  contains 

three  stacks  of  Ns  sensor  nodes  each,  we  consider  a  snapshot  of  the  complex  signal 

received  by  each  node  in  a  particular  stack,  which  we  model  as  the  true  signal  received  at 
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that  stack  with  phase  perturbations  due  to  position  errors.  By  using  the  mean  value  of  the 
in-phase  and  quadrature  components  of  all  the  measurements  from  a  particular  stack,  we 
can  estimate  the  true  signal  magnitude  and  phase  at  that  stack.  Since  position  errors  do 
not  affect  the  mean  value  of  the  phase  estimate,  this  method  will  recover  the  signal 
magnitude  and  phase  despite  the  presence  of  position  errors.  Given  this,  we  derive  an 
expression  for  the  mean  of  the  in-phase  and  quadrature  components  at  each  stack.  This 
expression  is  used  to  estimate  the  phase  and  magnitude  based  on  the  sample  mean  of  the 
received  complex  snapshots. 

Considering  an  ASLA  configured  network  with  uniform  position  errors,  we 
represent  a  snapshot  of  the  complex  signal  received  at  the  zth  node  of  the  center  stack  as 


Yc 


_  y  gifa  eiP{Sxisin0t  +SyiCOS0,  ) 


(114) 


where  i  =  l,...,Ns,  and  Vt  and  ^  are  the  signal’s  magnitude  and  phase,  respectively. 
Using  the  same  assumptions  used  to  derive  (106),  we  express  (1 14)  as 

r.j  =  U<*"vl  (ii5) 

where  p£ ,  =  J3SV ,  is  a  uniform  random  variable  representing  the  zth  node’s  phase 

perturbation  due  to  position  errors  in  the  v-direction.  The  resulting  complex  signal  is 
then  represented  as  two  random  variables  [85] 

rc,i=Yc,.i+jYcQ.i>  (H6) 

where  y  .  is  the  zth  node’s  in-phase  component  expressed  as 

Yc„i  =  K,cos  ($+/?,.,■)  (117) 

and  yC(> is  the  z'th  node’s  quadrature  component  expressed  as 
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YcQ,i  =  Ksin(fi+Psj)- 


(118) 


We  consider  ;/  and  7  as  two  random  variables  representing  a  given  stack’s  received 

in-phase  and  quadrature  components,  respectively.  The  random  variables  yC/ .  and  y  , 

can  be  interpreted  as  the  zlh  instance  of  the  true  and  /  that  are  perturbed  by  node 

position  errors.  Also,  the  mean  values  of  all  the  instances  would  then  be  the  sample 
means  ofyCi  and  y  .  Since  we  have  shown  that  position  errors  do  not  affect  the  mean 

value  of  the  signal  phase,  we  can  use  the  sample  mean  of  yC/  and  to  derive  an 
estimate  of  the  true  signal  phase.  Overall,  by  deriving  an  expression  for  the  mean  of  y 
and  y ,  we  can  also  determine  a  way  to  estimate  the  true  signal  phase. 

With  a  goal  of  derive  an  expression  for  the  mean  of  y  and  for  a  given  stack, 
we  begin  by  deriving  the  probability  density  function  of  y  .  By  substituting  ps  =  (j)t  +  ps 
in  (116),  we  get 


y  =Vpm{ps). 


(119) 


The  probability  density  function  of  ps  can  be  expressed  as 


•4(a) 


^  -  PSp  -  Ps  -<t>t+  PSp  • 


(120) 


We  next  substitute  ks  =  sin  ( p. )  in  (1 19)  to  get 


Ycq  =  V,Ks 


(121) 


The  probability  density  function  of  the  random  variable  ks  is  expressed  as 
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/.>.)= rr^-r 


(122) 


where 


g  1(rs)  =  5m^1/f^ 


(123) 


<i/r  / dp c  ,,  ,  =  cos(  sm  1  r,l 

5  s  A=g  (*>)  V  s/ 


=  a/1  —  C  at„ 


(124) 


for  sin  (V/)  - j  <  v  <  sin  (r/i  +  //()„ ) .  Finally,  the  probability  density  function  of  the 


random  variable  yc  =  is  expressed  as 


/re  (fcj- ^ /ks  f 


(125) 


=  iL  Itt-2  1  2  '  V‘sin  fa  -  ^ *  Ksin  (f,  +PSp). 


2^PJvt2-rCn 


We  then  account  for  measurement  noise  using  the  expression 


h  Yc,  +^AWGN-> 


(126) 


Qs=7c+Qa 


(127) 


where  I A  WGN  and  QAWGN  are  independent  Gaussian  random  variables,  representing  the 
complex  additive  white  Gaussian  noise.  Using  these  expressions,  we  derive  an  estimate 
of  V,  and  (f>t  from  the  mean  of  all  the  center  stack’s  complex  in-phase  and  quadrature 
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measurements.  The  mean  of  each  stack’s  in-phase  and  quadrature  components  can  be 
expressed  as 


Mc,=E{I,}  =  E{rCi}  +  E{IAWm),  (128) 

He,  =E{Q,}  =  E{rCa)  +  E{QAm„}.  (129) 


With  regard  to  the  quadrature  component,  since  QiWGN  is  zero  mean,  //r  can  be  re¬ 
written  as 


'B0  1 


rCn 


Saq  2 a, 


K-rl 


^d7Cn  + 


f- 

Jab  2an 


rCn 


-Bn  1 


Ycn 


2 a 


rdYcn> 


:dyc,. 


K-rl 


^-ap<  +  S 


(130) 


W>§  ~aP 


where  ap  =  (58 p ,  AQ  =  Vtsin  (<j>t  -  a p ) ,  and  B0  =  Vtsin  [<(  +  ap ) .  We  define  the  mean  of 
the  in-phase  component  similarly  as 
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Me,  E\yCi 


j 'Bi 

J  A/ 


rB, 

J  A/ 


1 

yCl 

lap . 

JA-A, 

1 

2ap  - 

h2-rl, 

1 

y Cj 

2ap . 

h2-A 

d/c,’ 


~aP 


,  rK  1 

dr^UTa 


rc, 


2a,fi, 


2  A 


dYCl> 


n 

- a„< 

2 


p  \Tt 


.  n 

<  —  +  aD 
2  p 


drCl> 


W>§ 


(131) 


where  At  =  F,cos(y,  -ap)  and  Bl  =  Vtcos (<pt  +ap)  ■  By  solving  (130)  for  (f)t ,  we  obtain 
an  expression  for  the  true  signal  phase  that  is  a  function  of  /j,c  ,  Vt  and  ap .  Then  using 
the  sample  mean  juc^  in  place  of  the  true  mean  jli(  ^  ,  we  obtain  the  phase  estimate  (j)t . 
The  derivation  of  the  phase  estimate  (j),  for  the  case  of  \<j)t  |  <  n  /  2  -  a p  is  shown  as 


Mc„  = 


[BQ  1  YCn 


0  jAe  2a, 


dyCn . 


2  Q 


2a. , 


(132) 


AQ 


2a, . 


By  squaring  both  sides  and  expanding  Aq  and  Bq  terms,  we  get 
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Expanding  the  right  hand  side  and  applying  the  trigonometric  identity 
cos2  0  -  sin2  0  =  cos  20  ,  we  get 


(2apVcg)  =E,2(l  +  cos(2  </>,)  cos  (; 2°p )) 

~2yjvt4cos2  (</>,  -  ap  )cos2  ($  +  ap ) 


(134) 


By  taking  V4  out  of  the  square  root  and  bringing  Vt2  to  the  left  hand  side,  we  get 


{2apMcQ )  =E,2(l  +  cos(2^)  cos  (2a«.)) 

-2E,2  |cos($  -ap)cos($  +  ap) 


(135) 


By  applying  trigonometric  identity  cos  x  cosy  =  (1/2)  [cos  (x  +  y)  +  cos  (x-y)]  to  the 
absolute  value  terms,  we  get 

V 

=  1  +  cos(2^,)cos[2a;j)-  cos(2ap)+  cos (2^)  .  (136) 

) 


We  perform  the  same  derivation  for  the  mean  of  the  in-phase  component  to  obtain 


=  l-cos(2^)cos(2arp  )-  cos(2cf;^-cos(2^,)  .  (137) 


Taking  the  ratio  of  (136)  and  (137)  yields 

K 

fic, 


1  +  cos  (2 (f)t )  cos  (2a p  )  -  cos  (2ap)+cos(2  $) 
1  -  cos  (2<f)t )  cos  [2a p )  -  cos  [2a p  )  -  cos  (2$ ) 


(138) 


'2 

l  K  ' 
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To  solve  for  ,  we  substitute,  jus  =  ( juCq  /  juC/  )  ,  a  =  cos(2ap)  and  b  =  cos(2^,)in  (138) 
to  get 


\  +  ab-\a  +  b\ 
\-ab-\a-b\ 


(139) 


By  taking  only  the  positive  solutions  for  the  absolute  value  terms  and  solving  for  b ,  we 
get 

b=]—^.  (140) 

l  +  /k 


Since  b  =  cos  (20, ) ,  we  can  solve  (140)  for  the  phase  estimate 

\ 

^=(l/2)cos_1  — —  .  (141) 

V  1  t*s  y 


For  the  other  two  cases  of  (130),  repeating  the  above  steps  yields  the  same  expression  for 
the  phase  estimate  <f)t .  By  substituting  (f)t  into  (137)  and  solving  for  Vt,  we  get  the  signal 
magnitude  estimate 
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With  each  stack  in  the  array  using  the  same  estimation  technique,  we  express  their 
complex  signal  estimates  as 

&n=K,,exP  (M,,)’  n  =  1,2,3,  (143) 

where  Vtn  and  (j)t  n  are  the  left,  right  and  center  stack’s  magnitude  and  phase  estimates, 
respectively.  Each  stack’s  Q„  is  then  treated  as  the  complex  output  of  a  single  node  in  a 
phase  shift  beamforming  array.  The  final  main  beam  response  is  then  expressed  as 

n,=Zn,.  (144) 

n= 1 

1,  Signal  Collection  in  the  Presence  of  Thermal  Noise  and  Uniform 
Sensor  Position  Error 

To  demonstrate  the  performance  of  our  signal  collection  and  estimation  scheme, 
we  present  the  results  of  five  Monte  Carlo  simulations  in  Figure  41  through  Figure  45. 
The  results  of  these  simulations  focus  on  the  root  mean-square  error  for  signal  magnitude 

(£)  and  phase  estimate.  For  all  five  simulations,  X  =  1  m,  the  source  emitter  is 

located  at  cot  =  [2000  m,  45  deg]T ,  r  =  200  meters,  and  all  data  points  are  the  result  of 
10,000  trials.  Note  a  key  parameter  in  the  following  simulations  is  signal  to  noise  ratio 
(SNR)  in  units  of  dB.  This  quantity  is  a  measure  of  the  ratio  between  signal  power  Pr 

and  noise  power  Pn  and  is  expressed  as  SNR  =  101og10  ( Pr  / Pn) . 

A  comparison  of  the  scheme’s  and  ^  as  a  function  of  per  node  SNR  at 

different  RDOA  noise  variance  o)U}  levels  are  shown  in  Figure  41  and  Figure  42, 

respectively.  The  target  bearing  estimate  6t  is  derived  using  (75).  From  the  results,  we 

see  that  the  effect  of  the  per-node  SNR  on  and  is  not  as  strong  as  that  of  the  6t 

estimate.  This  suggests  that  the  estimator  is  heavily  dependent  on  the  accuracy  of  the 
localization. 
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Figure  41.  Root  mean-square  error  at  different  levels  of  RDOA  noise  versus 
SNR  per  node.  Results  based  on  10,000  trials  with  Ns  =  10 ,  ra  =  200  m , 
8 p  =  .05  m,  and  a>t  =  [2000  m,  45  deg]T . 


— x —  RDOA  noise  (-1 5  dB) 

— B —  RDOA  noise  (-10  dB) 

—  RDOA  noise  (-5  dB) 

RDOA  noise  (0  dB) 

Figure  42.  Root  mean-square  error  ^  at  different  levels  of  RDOA  noise  versus 
SNR  per  node.  Results  based  on  10,000  trials  with  Ns  =  10 ,  r  =  200  m , 
8p  =  .05  m,  and  cot  =  [2000  m,  45  deg]T . 
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Two  plots  comparing  the  scheme’s  ^  as  a  function  of  maximum  unifonn 
position  error  Sp  and  the  number  of  nodes  per  stack  Ns  as  a  function  of  a]a)  are  shown  in 
Figure  43  and  Figure  44,  respectively.  In  Figure  43,  we  see  that  uniform  position  error 
doesn’t  affect  ^  until  values  of  ^p>0.1m  in  the  low  o\D  case =  -10  dB) . 
Conversely,  it  causes  a  near  linear  increase  in  ^  for  all  values  of  8  simulated.  In 
Figure  44,  the  results  show  that  the  benefits  of  increasing  Ns  diminishes  when 
Ns  >  10  nodes  .  Unlike  the  results  from  Figure  21,  where  increasing  the  number  of  nodes 

improved  performance,  here  increasing  the  number  of  nodes  per  stack  is  less  effective 
above  a  certain  value. 


Figure  43.  Root  mean-square  error  ^  at  different  levels  of  RDOA  noise  under 
increase  maxed  uniform  position  error  8p  .  Results  based  on  10,000  trials 
with,  Ns  =  10,  r  =  200  m ,  SNR  =  5  dB,  and  (ot  =  [2000  m,  45  deg]T . 


Figure  44.  Root  mean-square  error  ^  at  different  levels  of  RDOA  noise  under 
increase  maxed  uniform  position  error  Ns .  Results  based  on  10,000  trials 
with,  8p=\/\Qm,ra=  200  m ,  SNR  =  5  dB,  and  cot  =  [2000  m,  45  deg]T . 
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A  plot  of  the  proposed  scheme’s  E,v  as  a  function  of  per  node  SNR  with  and 
without  signal  estimation  is  shown  in  Figure  45.  Each  case  was  also  evaluated  at  two 
different  values  of  alD .  The  “with  estimation”  case  represents  our  proposed  scheme,  i.e., 

collaborative  beamforming  in  conjunction  with  the  proposed  signal  estimator.  The  “no¬ 
estimation”  case  represents  collaborative  beamforming  without  the  use  of  signal 
estimation.  In  this  simulation,  the  source  emitter  is  located  at®,  =[2000  m,  45  deg]T, 
Ns  - 10 ,  A  =  1  m ,  S  —  1  / 10  m ,  and  all  data  points  are  the  result  of  10,000  trials.  In  the 
o\D  =-15  dBcase,  we  see  that  the  benefits  of  the  signal  estimation  are  evident  at  all 
values  of  SNR  simulated.  As  evident  in  the  ct2rd  =  -5  dB  case  where  the  scheme 
underperforms  the  no-estimation  case,  we  again  see  that  the  scheme  is  highly  dependent 
on  the  accurarcy  of  the  localization  phase.  This  dependency  is  derived  from  the  scheme’s 
assumption  that  the  target  signal’s  bearing  is  perfectly  normal  to  the  array,  i.e.,  the 
6 t  =  0  deg,  when  in  actuallity  6t  it  is  a  small  random  quantity  based  on  the  accuracy  of 
the  localization. 


Figure  45.  Signal  of  proposed  scheme  compared  with  the  no-estimation  case 
versus  SNR  per  node  values.  Results  based  on  10,000  trials  with 
=  1  / 10  m ,  ra  =  200  m ,  and  a>t  =  [2000  m,  45  deg]T . 
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An  illustration  of  the  mean  value  of  the  normalized  gain  of  the  array  as  a  function 
of  signal  frequency  for  different  values  of  Sp  is  shown  in  Figure  46.  For  comparison,  we 
also  include  the  no-estimation  case.  In  this  simulation,  the  source  emitter  is  located  at 
cot  =  [2000  m,  45  deg]T ,  Ns  =  10 ,  A  =  1  m ,  Sp  =  1  / 10  m ,  o ^  =  0  dB,  and  all  data  points 
are  the  result  of  10,000  trials.  As  frequency  increases,  so  does  J3  .  From  (115),  we  see 
that  phase  error  due  to  Sp  is  proporitional  to  ft .  As  (3  =  2k  fc  /  c ,  an  increase  in 

frequency  will  magnify  the  negative  effects  of  the  position  errors.  From  the  results,  we 
see  that  allowing  these  position  errors  into  beamforming  (no-estimation  case)  results  in  a 
response  similar  to  that  in  Figure  36  where  the  main  beam  gain  is  significantly  reduced  as 
position  errors  increased.  When  we  incorporate  the  proposed  signal  estimator,  an 
improvement  over  the  no-estimation  case  is  apparent  for  all  simulated  frequency  values. 
This  suggests  that  the  inclusion  of  signal  estimation  in  collaborative  beamforming  can 
effectively  mitigate  array  phase  pertabations  due  to  the  uniformly  distributed  positioin 
errors. 


Figure  46.  Average  signal  normalized  gain  at  different  levels  of  RDOA  noise 
versus  signal  frequencies.  Results  based  on  10,000  trials  with  cr^  =  0  dB, 

8p  =  1  / 10  m ,  ra  =  200  m ,  and  (ot  =  [2000  m,  45  deg]T . 
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D.  SIGNAL  COLLECTION  FROM  AIRBORNE  SYMMETRIC  LINE  ARRAY 
NETWORK  IN  THE  PRESENCE  OF  GAUSSIAN  POSITION  ERROR 

In  the  previous  section,  we  developed  a  signal  estimator  to  combat  the  effects  of 
uniformly  distributed  position  errors.  In  this  section,  we  derive  a  similar  estimator  for  the 
case  of  Gaussian  distributed  position  errors. 

In  a  fashion  similar  to  (116),  a  snapshot  of  the  complex  signal  received  at  the  7th 
node  of  the  center  stack  is  represented  as 

nc,i=nc„i+ncg,i>  O45) 

where 

nc„i  =C(cos  ($+/?„),  (146) 

with  pn  -  P5y  with  S  being  zero  mean  Gaussian  position  error  with  a  variance  of  cr 
(see  Chapter  III)  and 

ncQ,i  =  JKs'n  (</>,+  P,,)-  (147) 

Similar  to  the  uniform  position  error  case,  we  will  derive  an  expression  for  the  mean 
value  of  the  in-phase  and  quadrature  components  for  each  stack  and  then  use  these 
expressions  to  derive  an  estimate  of  the  signal  phase.  Also  similar  to  the  uniform  random 
variable  case,  each  node’s  nCi .  and  nc^ ;  in  a  stack  can  be  interpreted  as  the  7th  instance  of 

the  in-phase  random  variable  nc  and  the  quadrature  random  variable  nc^  . 

Focusing  on  the  in-phase  component  first,  the  expression  (f>t  +  pn  can  be  expressed 
as  a  Gaussian  random  variable  8  with  a  mean  of 

Pg=<f>,  (148) 

and  variance  of 
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(149) 


By  substituting  8  into  nc  ,  we  get 


(150) 


To  derive  the  resulting  probability  density  function  of  nC) ,  we  first  define  its  cumulative 
distribution  function  as 


(nQ)  =  Pr[g(^)<nC/]=  J  fUg(Sg)dSg 


(151) 


g(^)inc 


where  g  ( 8 g  j  =  cos  ( 8 g  j .  Since  a  Gaussian  random  variable  is  define  between  -oo  and  oo, 
5a  can  be  defined  by  the  inverse  transform  of  g  ( 8g  j  over  the  interval  [84] 


2 nk  +  cos  1  [nCi  )<8g  <  2x(l  +  k )  -  cos  1  [nCi  j . 


(152) 


The  probability  density  function  of  is  then  expressed  as  [81],  [86] 


/»C,("C,)  =  E 


2  nk  +  cos 


1  (nc;  ))  +  /„  (2^(1  +  k)-  cos1  [nCj ))) ,  (153) 


where 


dnCi 


(154) 
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By  expanding  (153),  we  get  the  final  expression  for  the  probability  density  function  of 


nCj  as 


oo 

K  {nc,)=  X 


1 

1 

~  llCi 

1 

exp 


A -6 


2  a 


+  exp 


s  y 


Br~*t 


W 


2(7 


(155) 


S  2) 


where  Ay  =  2n k  +  cos  1  (nC; )  and  By  =  2x(\  +  k)~  cos  1  [nCi )  . 

To  support  the  validity  of  this  expression,  a  comparison  of  the  histogram  of 
n()  =  f) cos  (<:>'„  j  with  the  response  of  (155)  over  the  same  interval  is  shown  in  Figure  47. 

For  this  simulation,  10,000  instances  of  Su  generated  to  form  a  histogram  with  50  bins, 
cr  =  {2n  /  5  )2 ,  and  fNr  [nCi  j  is  calculated  with  £  =  0,1.  As  observed  from  the  results,  the 
simulated  histogram  closely  follows  the  theoretical  values  of  (155). 


Figure  47.  Comparison  of  measured  and  theoretical  probability  density  function 
(pdf)  of  nCi .  Signal  phase  =  0  and  measured  pdf  based  on  histogram 
using  10,000  points  and  50  bins. 
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1.  Mean  Value  of  the  In-Phase  and  Quadrature  Components  with 
Gaussian  Position  Errors 

With  the  probability  density  function  of  nCj  defined  as  a  function  of  (j)t  and  cr2 , 
we  can  derive  an  estimate  of  4>l  at  a  given  stack  based  on  the  sample  mean  of  n(/  and 
knowledge  of  cr2 .  Since  it  is  difficult  to  derive  an  expression  for  the  mean  using  the 
probability  density  function  of  nCi ,  we  take  an  alternative  approach  using  Euler’s  fonnula 

e'dg  =  cos 8g  +  ysin^g  .  (156) 

We  will  first  consider  the  case  where  8g  is  Gaussian  with  a  mean  value  of  zero  and  a 
variance  of  cr2 .  Note  that  the  mean  value  of  8  u  is  indeed  equal  to  (f>t  (see  (148)),  but  we 

will  first  examine  this  simpler  zero  mean  case  in  order  to  derive  the  non-zero  mean  case. 
The  mean  of  the  cosine  term  is  then  given  as  [87] 

E  jcosci,}  =  exp(-cr2  /  2), 
and 

E  jsint>g  j  =  0. 

Using  the  trigonometric  identity 

cos2^=|(l  +  cos(2^)),  (159) 

the  variance  of  cos  8a  can  be  expressed  as 

var|cos^g|  =E'|cos2^gj-E'|cos^g|  .  (160) 


(157) 


(158) 
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By  substituting  (157)  and  (159)  into  (160),  we  get 


var 


{C0S<^}  =  T 


f  /o  _  \2  /o  A  _2  I 


,-K)/2 


1  +  e 

\  j 


-  —\\-e 


(161) 


Similarly,  the  variance  of  sin£  is  given  by 


var | sin <5^ |  =  fsjsin2f>g}-fs{sin<7gj  =  —  ll-e  2°s 


(162) 


Finally  to  consider  the  case  where  Sg  is  not  zero  mean,  where  Sg  is  a  Gaussian  random 
variable  with  mean  //,,  =  (f)t  and  variance  equal  to  cr  .  If  this  is  true,  then  S0  -  Sg  -jug 
where  S0  is  the  zero  mean  case  of  §a .  Thus,  using  the  trigonometric  identity 
cos(*  +  y)  =  cos  v  cos  y  -  sin  *  sin  y  ,  the  mean  of  cos(t)g  )  can  be  expressed  as 

£■  {cost?,}  =  £ jcos(c)0  +  /us  jj  =  Zr|cosf)0cos//,,  -sin£0sin//g}  =  e_<Tg/2cos /ug.  (163) 
Similarly, 

fi  {sin  j  =  £  jsin(e>0  +  =  e_c7®/2sin  jug.  (164) 

Since  nC]  =  cos(<7„  j ,  the  mean  of  nC;  is  given  as 

E{nc  }=e  ag/2cos/ug.  (165) 

Following  the  same  derivation,  the  mean  of  nc  is  expressed  as 

E{ncQ}  =  e^l2s[nJug-  (166) 
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2.  Phase  Estimates  Using  Sampled  Base  Band  Signals 

Given  that  the  signal  phase  is  perturbed  by  sensor  position  errors,  as  shown  in 
Chapter  V  section  B,  the  mean  value  of  the  signal’s  phase  remains  largely  unaffected. 
This  makes  the  sample  mean  of  a  given  stack’s  signal  phase  a  suitable  parameter  from 
which  to  estimate  the  true  signal  phase. 

One  means  of  estimating  the  signal  phase  is  by  solving  (163)  or  (165)  for  (j)n  to 
obtain  the  phase  estimate  for  the  center  stack 

kc  =  cos_1  (A,q  / e~as/2j  =  she1  (//„Ce  / e"V2) ,  (167) 


where  jun<  and  jiu<  are  the  signal’s  in-phase  and  quadrature  component’s  sample  mean 

for  the  center  stack.  This  method,  although  valid,  only  uses  the  samples  from  either 
component.  The  preferred  approach  is  to  use  both  base  band  signal  samples  by  taking  the 
ratio  of  (166)  and  (165)  as 


E{ncQ]  _  e~a,'2sm Ms,  _  sitVC,; 
EM  e'^/2cos ns  cos Hst  ' 


Since  ju5  =  (/)!  c ,  an  estimate  of  <pt  c  is  expressed  as 


kc  = tan 


fju  > 


v  ^  j 


(168) 


(169) 


This  same  phase  estimation  method  is  then  used  for  the  other  two  stacks  in  the  array.  The 
final  signal  phase  estimate^,  is  then  the  average  of  these  three  estimates.  Since  the 
ASLA  array  is  symmetrical  about  the  y-axis,  taking  the  average  of  the  estimates  will  also 


98 


compensate  for  the  minor  errors  in  the  arrays  reorientation,  i.e.,  errors  in  the  estimate  6t . 
The  estimated  complex  signal  for  the  nth  stack  is  then  expressed  as 

=exp(.M),  (170) 

We  then  sum  the  estimates  of  all  three  stacks  to  obtain  the  final  main  beam  response  as 

4 =Z4,  (i7i) 

n= 1 


3.  Performance  of  the  Signal  Collection  Scheme  in  the  Presence  of 
Gaussian  Sensor  Position  Errors 

To  evaluate  the  proposed  scheme’s  performance  in  the  presence  of  Gaussian 
sensor  position  errors,  we  conducted  a  series  of  Monte  Carlo  simulations.  The  results  of 
these  simulations  focus  on  the  root  mean-square  error  for  signal  magnitude  (<^v)  and 

phase  estimate.  For  comparison,  the  no-estimation  case,  i.e.,  beamforming  without 
signal  estimation,  is  also  shown.  For  each  of  these  simulations,  A  =  1  m ,  the  source 
emitter  is  located  at  cu,  =  [2000  m,  15  deg]T ,  r  =  200  meters,  and  Ns  =10.  All  data 
points  are  the  result  of  10,000  trials. 

A  comparison  of  the  scheme’s  root  mean-square  error  of  normalized  magnitude 
t,v  with  and  without  signal  estimation  as  a  function  of  the  standard  deviation  of  Gaussian 

position  error  <Jg  is  shown  in  Figure  48.  In  this  simulation,  to  isolate  the  effect  position 

error,  the  target  signal’s  6t  is  known  a  priori.  From  the  results,  we  can  see  that  the 

estimated  case  outperforms  the  no-estimation  case  at  all  values  of  cr,  simulated.  This 

suggests  that  signal  estimation  provides  an  effective  means  of  reducing  the  effects  of 
Gaussian  position  errors. 
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Figure  48.  Root  mean-square  error  of  normalize  magnitude  versus  the  standard 
deviation  of  Gaussian  position  error  <r„ .  Results  based  on  10,000  trials  with 

Ns  =  10,  r  =  200  m ,  and  a>t  =  [2000  m,  15  deg]T . 

The  scheme’s  ^  as  a  function  of  Gaussian  positon  error  <Jg  is  shown  in  Figure 

49.  Similar  to  the  previous  simulation,  we  assumed  the  target  signal’s  6t  is  known  a 

priori.  From  the  results,  we  observe  no  noticeable  difference  in  the  ^  between  the 

proposed  signal  estimation  case  and  the  no-estimation  case.  Similar  to  the  uniform 
position  error  case,  this  suggests  that  the  proposed  signal  estimation  technique  does  not 
improve  the  accuracy  of  the  main  beam  signal  phase,  but  as  seen  in  Figure  48  does 
improve  its  magnitude  gain. 


100 


Figure  49.  Root  mean-square  error  phase  versus  the  standard  deviation  of 
Gaussian  position  error  <r„ .  Results  based  on  10,000  trials  with  Ns  =  10, 
r  =  200  m ,  and  cot  =  [2000  m,  15  deg]T . 

The  scheme’s  as  a  function  of  RDOA  noise  variance  (T2tD  =  c2<j2r  is  shown  in 
Figure  50.  In  this  simulation,  the  scheme  is  evaluated  at  two  levels  of  position  error  with 
cr„  =0  and  a,,  =1/5  m.  From  the  results,  we  can  see  that  the  proposed  scheme  and  no¬ 
estimation  case  have  the  same  performance  when  no  position  errors  are  present.  When 
position  errors  are  introduced,  we  see  that  the  proposed  scheme  outperform  the  no¬ 
estimation  case.  Of  interest,  we  see  that  the  scheme  also  outperforms  the  same  case  with 
no  position  error  for  values  of  RDOA  noise  greater  than  2  dB.  This  would  suggest  that  at 

a  given  a ^  level,  position  errors  become  less  of  a  factor  than  errors  in  6t . 
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RDOA  Noise  dB  (c 2ol) 

Figure  50.  Root  mean-square  error  of  normalize  magnitude  versus  RDOA  noise. 
Results  based  on  10,000  trials  with,  Ns  =  10,  ra  =  200  m ,  and 

co,  =  [2000  m,  15  degf. 


A  comparison  of  the  scheme’s  ^  as  a  function  of  a ^  is  shown  in  Figure  51. 
Similar  to  the  previous  simulation,  the  performance  was  evaluated  at  two  levels  of 
position  error,  aa  =  0  and  ag  =  1/5  m.  Similar  to  Figure  49,  we  can  see  that  the  proposed 

scheme  and  no-estimation  case  have  the  same  performance  with  or  without  position 
errors.  Also,  similar  to  the  previous  simulation,  we  see  that  the  cases  with  position  error 
beat  the  no  position  error  cases  at  values  of  RDOA  noise  greater  than  2  dB. 
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Figure  5 1 .  Root  mean-square  error  Phase  ^  versus  RDOA  noise.  Results  based  on 
10,000  trials  with,  Ns  =  10,  ra  =  200  m ,  and  cot  =  [2000  m,  15  deg]T . 

E.  GRATING  LOBE  SUPPRESSION  VIA  VIRTUAL  FILLING 

From  (55)  and  (56),  we  know  that  increasing  the  distance  r  between  the  outer 

and  center  stacks  (see  Figure  6)  will  provide  increased  localization  performance,  but 
increasing  this  distance  essentially  increases  the  inter-node  distance  d  in  (28).  If  the 
inter-node  distance  of  an  array  is  greater  than  X  /  2  ,  multiple  grating  lobes  will  appear  in 
the  array  factor  [71].  Since  grating  lobes  have  the  same  magnitude  response  as  the  main 
beam,  they  will  significantly  reduce  the  array’s  interference  rejection  capabilities. 

To  combat  the  effects  of  grating  lobes,  we  propose  the  use  of  virtual  filling. 
Given  that  the  directions  of  arrival  of  all  signals  are  known,  virtual  filling  interpolates  the 
response  of  a  virtual  uniform  linear  array  based  on  signal  estimates  from  the  real  array 
[76].  In  virtual  filling,  any  gap  between  array  nodes  larger  than  X  /  2  is  filled  in  with 
virtual  complex  signal  data. 
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With  the  ASLA  symmetrically  populated  only  along  the  x-axis  and  centered  on 
the  origin,  the  position  of  the  7?th  virtual  node  is  expressed  as 

=  ~ra  +  nvda  >  nv  =  0, 1, . .  .,iNTv  1  (172) 

where  the  virtual  inter-node  distance  da  =  A  /  2  and  Nv  is  the  total  number  of  virtual 
elements  given  as 

K=^+ 1.  (173) 

d„ 


To  fill  in  the  complex  data  of  the  virtual  array,  we  begin  by  populating  the  vector 
An  from  (51)  with  the  signal  estimates  Q;i.  Using  the  least-squares  method  described  in 
Chapter  II,  an  estimate  of  each  signal’s  magnitude  and  phase  is  obtained.  The  response 
of  the  nth  virtual  node  is  then  expressed  as 


^n  =  YVtje^eiP^9'* 


(174) 


where  NR  is  the  number  of  signals  present,  i  —  1, . . . ,  NR ,  and  Vti  ,</>tl,  0t ,  are  the  7th 

signal’s  magnitude,  phase,  and  bearing,  respectively.  Once  fdled  out,  the  ASLA  then 
becomes  a  virtual  uniform  linear  array.  Since  the  intemode  distance  is  now</l/2,  the 
presence  of  grating  lobes  is  significantly  reduced.  Furthermore,  various  tapering  methods 
(see  Chapter  II)  can  be  applied  to  manage  the  remaining  sidelobes. 

With  no  errors  present,  the  array  factor  of  an  Ns=5  ASLA  configured  network  is 
shown  in  Figure  52.  For  this  simulation,  A  =  1  m  and  ra  =  200  A  .  From  the  plot,  we  can 
see  that  with  ra  »  A  /  2 ,  there  are  numerous  grating  lobes  and  side  lobes  in  the  array 
factor. 
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Figure  52.  Normalized  Array  factor  of  an  Ns=5  ASLA  phase  shift  beamformer 
before  virtual  filling  process.  Simulation  based  on  ra  =  200  m ,  and 
cot  =  [2000  m,  45  deg]T . 

Now,  with  virtual  filling  and  tapering  applied,  the  array  factor  of  the  ASLA  is 
shown  in  Figure  53.  With  the  virtually  filled  array  having  no  inter-node  spacing  greater 
than  XI 2,  we  can  see  from  the  plot  that  all  the  grating  lobes  have  been  eliminated.  By 
applying  a  binomial  taper  (see  Chapter  II)  the  sidelobes  have  also  been  suppressed, 
resulting  in  a  strong  main  beam  lobe. 


Figure  53.  Normalized  Array  factor  of  an  Ns  =  5  ASLA  phase  shift  beamformer 
with  virtual  filling  process  with  binomial  tapering.  Simulation  based  on, 
r  =  200  m ,  and  a>t  =  [2000  m,  45  deg]T . 
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A  comparison  of  the  signal  estimation  scheme  with  and  without  virtual  filling  is 
shown  in  Figure  54  and  Figure  55.  In  this  simulation,  there  is  an  interfering  signal  in 
addition  to  the  target  signal,  A  =  1  m ,  uniform  position  error  Sp  =  1  /  20  in  ,  and  the  ASLA 

formation  consists  of  Ns  = 10  with  ra  =100  m.  The  target  signal’s  Gt—  0  deg,  Vs  =1, 
and  <f>s  =  60  deg .  The  interfering  signal’s  9t  -  55  deg ,  Vs  =  1 ,  and  </>s  —25  deg . 

The  magnitude  and  phase  estimate  values  of  the  main  beam  for  100  trials  without 
virtual  filling  are  shown  in  Figure  54.  From  the  results,  we  can  see  that  the  interfering 
signal  corrupts  the  target  signal,  thus  resulting  in  an  erroneous  main  beam  response.  In 
contrast,  the  main  beam  signal  values  using  virtual  filling  are  shown  in  Figure  55.  From 
the  results,  it  is  clear  that  the  virtual  filling  has  reduced  the  effects  of  the  interfering 
signal  and  has  improved  the  signal  estimate. 
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b)  Phase  estimate  without  virtual  filling 


Figure  54.  Signal  Estimation  without  virtual  filling.  Simulation  based  on, 
r  =  1 00  in  ,  target  signal’s  6t-0  deg ,  Vs  =  1 ,  and  <ps  =  60  deg .  The 

interfering  signal’s  0t  —  55  deg ,  Vs=l,  and  (f>s  =  25  deg .  The  true  signal 
phase  and  magnitude  is  shown  in  red. 
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a)  Magnitude  estimate  with  virtual  filling 


Figure  55.  Signal  estimation  with  virtual  filling.  Simulation  based  on,  r  =  100  m , 
target  signal’s  6t=  0  deg  ,  Vs  —  1 ,  and  <f>s  =  60  deg  .  The  interfering  signal’s 
6t  =  55  deg ,  Vs  =  1 ,  and  </>s  =  25  deg .  The  true  signal  phase  and  magnitude 

is  shown  in  red. 


A  comparison  of  the  main  beam’s  angular  root  mean  squared  error  ^  for  the 
proposed  scheme  with  and  without  virtual  filling  as  a  function  of  maximum  uniform 
position  error  8  is  shown  in  Figure  56.  In  this  simulation,  there  is  one  signal-of-interest 

and  one  interfering  signal,  the  target  signal’s  9t-0  deg,  Vt-l,  and  <f)t  =60  deg,  the 

interfering  signal’s  0t  =  55  deg ,  Vt  —  1 ,  <f)t  =  25  deg ,  and  ra  =  1 00  m  ;  and  all  data  points 

are  the  result  of  1,000  trials.  From  the  results,  we  see  in  the  no  virtual  filling  case  that  the 
benefits  of  the  signal  estimation  are  nullified  by  the  presence  of  an  interfering  signal.  In 
contrast,  when  virtual  filling  is  applied,  the  interfering  signal  is  significantly  reduced, 
thus  improving  the  signal  phase  estimate.  Overall,  this  result  demonstrates  the  proposed 
scheme’s  effectiveness  in  the  presence  of  both  position  errors  and  an  interfering  signal. 
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Figure  56.  Proposed  scheme’s  ^  with  virtual  filling  (blue)  and  without  (black) 
versus  maximum  uniform  position  error  5  in  the  presence  of  one 
interfering  signal.  In  this  simulation,  ra  =  100  m ,  the  target  signal’s 
6t  —  0  deg ,  Vt  =  1 ,  and  <f>t  =  60  deg ,  the  interfering  signal’s  6t=  55  deg , 

Vt=l,  and  <j)t  =  25  deg ,  and  all  data  points  are  the  result  of  1,000  trials. 

In  this  chapter,  we  studied  collaborative  signal  collection  in  the  presence  of  noise, 
sensor  position  errors,  and  interfering  signals.  We  developed  a  signal  estimation  method 
that  utilizes  sampled  array  data  and  knowledge  of  the  positional  error  statistics  to  combat 
these  errors.  We  further  enhanced  its  robustness  against  interfering  signals  by  using  a 
virtual  filling  and  tapering  techniques  to  minimize  the  array’s  grating  lobes  and  side  lobe 
response.  Using  simulation,  we  were  able  to  validate  the  proposed  scheme’s 
effectiveness  in  the  presence  of  such  errors  and  achieve  the  intended  objective  of  this 
dissertation. 
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VI.  CONCLUSION 


Performing  collaborative  beamforming  from  an  elevated,  mobile  wireless  sensor 
network  requires  a  coordinated  interplay  of  many  different  signal  processes  and 
technologies.  The  objective  of  this  research  was  to  maximize  signal  collection 
performance  in  the  presence  of  various  signal  and  sensor  related  errors  from  an  elevated, 
mobile  wireless  sensor  network.  To  accomplish  this  objective,  we  proposed  a  signal 
collection  scheme  that  exploits  an  elevated,  mobile  network  to  maximize  the 
collaborative  collection  of  a  target  signal. 

The  proposed  scheme  begins  with  an  accurate  source  localization  technique, 
which  is  used  to  detennine  the  signal’s  position  and  bearing.  This  infonnation  was  then 
used  by  a  collaborative  beamformer  to  maximize  the  collection  of  that  target  signal.  Both 
of  these  phases  are  enhanced  by  the  versatile  nature  of  an  elevated,  mobile  network.  The 
mobility  allows  for  the  reconfiguration  of  the  sensor  network’s  topology  to  help  create  an 
ideal  sensor-target  geometry.  This  geometry  allows  for  optimal  localization  estimates  to 
be  obtained.  The  ability  of  an  elevated  network  capable  of  creating  sensor  stacks  allows 
for  robust  signal  collection  operations  in  the  presence  of  uniform  and  Gaussian  sensor 
position  errors. 

In  the  localization  technique,  to  obtain  precise  localization  in  the  presence  of 
noise  and  uniform  position  errors,  we  proposed  the  use  of  two  sequential  location 
estimators.  This  technique  consists  of  an  initial  weighted  least-squares  estimate  followed 
by  a  maximum-likelihood  estimate.  The  second  estimate  uses  the  initial  estimate  to 
reorient  the  network.  This  reorientation  creates  an  ideal  sensor  to  target  geometry  [17]. 
This  new  orientation  minimizes  the  geometric  dilution  of  precision  [52],  minimizing  the 
maximum-likelihood  estimator’s  error  variance.  Simulation  results  showed  the  proposed 
localization  technique  to  be  efficient,  i.e.,  the  variance  of  the  estimation  error  approaches 
the  Cramer-Rao  lower  bound. 

To  enhance  the  localization  technique’s  performance,  we  developed  a 
measurement  outlier  rejection  process  that  mitigates  the  effects  of  measurement  and 
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sensor  position  errors.  This  technique  uses  a  combination  of  single  case  diagnostics  and 
the  Mahalanobis  distance  to  identify  and  remove  erroneous  measurements  from  the  least- 
squares  estimate.  Through  simulation,  we  demonstrated  this  technique  to  be  effective  in 
the  presence  of  both  measurement  and  position  errors.  For  the  case  where  only  position 
errors  are  present,  the  technique  was  found  to  be  effective  only  for  high  levels  of  position 
errors. 

For  the  signal  collection  technique,  we  began  by  analyzing  the  effects  of  sensor 
position  errors  on  the  array  factor  and  its  main  beam  gain.  We  then  derived  an 
expression  for  the  effects  of  both  unifonn  and  Gaussian  position  errors  on  the  array’s 
main  beam  gain.  Through  simulation,  we  validated  our  results  and  showed  that  the  mean 
value  of  the  main  beam  signal  phase  was  unaffected  by  position  errors.  To  enhance  the 
signal  collection  perfonnance  in  the  presents  of  these  errors,  we  developed  a  signal 
estimator  for  both  unifonn  and  Gaussian  position  errors.  Simulation  results  yielded  an 
improvement  in  array  gain  of  approximately  37  percent  for  the  standard  deviation  of 
position  error  values  greater  than  0.4  m. 

Taking  advantage  of  the  concept  of  a  unique  elevated,  mobile  wireless  sensor 
network  realized  through  the  use  of  multirotor  UAVs,  our  scheme  used  two  existing 
localization  techniques  to  deliver  a  precise  source  location  estimate.  Using  this 
information  with  statistical  knowledge  of  the  sensor  position  errors  in  a  signal  estimator, 
we  derived  a  novel  collaborative  signal  collection  scheme.  This  scheme  was  shown  to  be 
capable  of  collecting  and  amplifying  a  target  signal  in  the  presence  of  such  errors.  With 
all  these  techniques  in  concert,  the  objective  of  this  dissertation  to  maximize  signal 
collection  in  the  presence  of  various  signal-  and  sensor-related  errors  was  achieved. 

A.  SIGNIFICANT  CONTRIBUTIONS 

With  the  objective  to  maximize  signal  collection  performance  in  the  presence  of 
various  signal  and  sensor  related  errors,  a  novel  scheme  for  robust  signal  collection  from 
an  elevated,  mobile  network  was  proposed. 

To  enable  the  signal  collection,  a  localization  step  was  required.  To  accomplish 
this  task,  a  two-stage  localization  technique  was  proposed.  Using  an  optimal  ASLA 
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sensor  formation  for  target  localization,  we  minimized  the  network’s  geometric  dilution 
of  precision,  effectively  lowering  the  estimate’s  Cramer-Rao  lower  bound.  By  utilizing 
two  unique  location  estimators  in  sequence,  we  created  a  hyperbolic  location  estimate 
capable  of  approaching  the  Cramer-Rao  lower  bound. 

To  enhance  the  localization  robustness,  a  measurement  outlier  rejection  process 
was  proposed.  Using  the  Mahalanobis  distance  as  a  measure  of  solution  divergence,  we 
were  able  to  identify  and  remove  erroneous  measurements  from  the  location  estimate.  As 
a  result,  the  localization  estimation  technique  demonstrated  increased  resilience  to  both 
measurement  and  sensor  position  errors. 

To  better  understand  the  effects  of  position  errors  in  collaborative  beamfonning, 
an  expression  for  uniform  and  Gaussian  sensor  position  error  effects  on  an  ASLA 
formation’s  main  beam  gain  were  developed.  The  main  result  was  the  validation  that 
position  errors  in  the  ASLA  formation  have  no  effect  on  the  mean  value  of  the  signal 
phase,  which  was  instrumental  in  the  development  of  a  novel  signal  estimator  used  to 
compensate  for  such  errors. 

Finally,  a  signal  collection  scheme  that  combines  a  novel  signal  estimator  and 
collaborative  beamforming  to  successfully  collect  a  target  signal  despite  the  effects  of 
sensor  position  errors  was  proposed.  Compensating  for  both  unifonn  and  Gaussian 
errors,  the  signal  estimator  uses  sampled  array  data,  sensor  stack  formations,  and 
knowledge  of  the  positional  error  statistics  to  accomplish  this  task.  Through  simulation 
we  were  able  to  demonstrated  significant  improvement  over  the  conventional 
beamforming.  To  further  enhance  the  collection  capability  for  noise  and  interference 
rejection,  we  also  adapted  the  use  of  virtual  filling  and  tapering  to  manage  the  array’s 
grating  lobes  and  sidelobes. 

B.  FUTURE  RESEARCH 

With  the  research  and  development  of  localization  and  signal  collection 
consistently  changing  and  growing,  there  are  a  number  of  possible  extensions  to  the 
research  presented  in  this  dissertation. 
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The  key  error  parameter  evaluated  in  this  research  was  sensor  position  errors 
modeled  as  a  uniform  or  Gaussian  random  variable.  The  assumption  of  specific 
probability  density  functions  affects  many  aspects  of  each  of  the  techniques  and  methods 
developed  in  this  research.  Further  research  focused  on  a  general  approach  to  treat 
position  errors  independent  of  a  probability  distribution  should  be  considered  in  order  to 
fully  investigate  their  effects  on  signal  collection. 

The  robust  source  localization  technique  presented  here  focused  on  a  single  case 
diagnostic.  This  method  assumes  that  single  data  points  act  alone  to  offset  the  true 
solution.  It  has  been  shown  that  multiple  points  in  unison  may  also  have  a  stronger  effect 
than  an  individual  data  point  [18].  Given  this,  exploring  the  use  of  outlier  subset 
rejection  in  the  context  of  hyperbolic  localization  may  lead  to  more  beneficial  results. 

Both  source  localization  and  collection  operations  presented  here  assume  a  line- 
of-sight  signal  path  between  each  sensor  and  the  transmitter.  Further  research  is  needed 
to  adapt  the  scheme  for  use  in  a  non-line-of-sight  enviromnent  where  there  is  also  a 
multipath  component  to  each  signal  path. 

With  all  wireless  signal  processing  operations,  there  must  be  a  consideration  for 
the  wireless  communication  scheme  as  well.  Further  research  is  required  to  evaluate 
signal  collection  from  an  elevated,  mobile  network  in  the  context  of  standard  digital 
wireless  modulation  schemes,  such  as  binary  phase-shift  keying,  quadrature  phase-shift 
keying,  and  quadrature  amplitude  modulation.  Additional  effort  should  be  devoted  to  the 
orthogonal  frequency-division  multiple  access  signals  [85]. 
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APPENDIX.  SELECTED  SIMULATION  CODE 


Script  Name:  Two-Stage  Localization  Versus  RDOA  Noise 

clc;  close  all;  clear  all;warning  off;  %  Program 
initialization . 

L=10000;  %  Number  of 

ensemble  runs . 

c  =  3e8;  %signal  propagation  speed; 

LL  =  400;  %deployment  diameter 

r  =  2500;  %range  from  origin  to  source 

theta_tx  =  degtorad ( 60 ) ;  %This  is  from  the  x-axis 
uo  =  [r*cos (theta_tx)  r*sin ( theta_tx) ] ' ; 

%  True  source  position. 

NN=5;  %Number  of  nodes  per  stacks,  total  sensor  =  3*NN 
x— [  0*ones(l,NN)  - (LL/2 ) *ones ( 1 , NN)  (LL/2 ) *ones ( 1 , NN) ] ; 

%  True  sensor  position  matrix. 
y=zeros ( 1 , length (x) ) ; 

S= [x;  y]  ; 

M  =  size  (S, 2);  %  Number  of  sensors. 

LM=length (x) ; 

N  =  size  (S, 1) ; 
localization . 

MLE_runs  =  1; 

ro  =  sqrt ( sum ( (uo*ones ( 1 , M) -S) . A2 ) ) ' ; 
sensor  ranges 
rdo  =  ro (2 : end) -ro ( 1 ) ; 

R  =  (eye (M-l ) Tones (M-l )) /2 ;  %  covariance 

structure  of  TDOA 

R2  =  (eye (LM-1 ) Tones (LM-1 )) /2 ; 

NsePwrVecdB=-10 : 1 : 5; 

TSweighted  least-squares=  zeros  (2,1); 

Linear_MLE  =  zeros  (2,1); 

fprintf (' Simulation  in  progress'); 

for  nseldx=l : length (NsePwrVecdB) ,  %  vary  measurement 

noise  level 

fprintf  ('.'); 

nsePwr  =  1 0 A (NsePwrVecdB (nseldx) / 1 0 ) ; 

Q  =  nsePwr  *  R;  %  Covariance  matrix 

of  TDOA  (range  difference)  noise 
Q2  =  nsePwr  *  R2 ; 

crib (nseldx) =sqrt (trace (TDOALocCRLB^lin (S, uo, Q) ) ) ; 
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%  Dimension  of 


%  True  source- 


Linear_range  = 

(  (  (1/nsePwr) *cA2* (1/ (c* (r) A2) ) A2* (LL/2) A4  * (LM/8) ) A (-1) ) ; 
crlb_LM (nseldx) =sqrt (Linear_range) ; 

o, 

o 

Linear_bearing  =  ( ( (2/nsePwr) *cA2* (1/ (cA2) ) * (LL/2) A2 
* (LM) ) A (-1)) ; 

crlb_LM_bearing (nseldx) =sqrt (Linear_bearing) ; 

SimulationMSEl  =  0; 

SimulationMSE2  =  0; 

SimulationMSE3  =  0; 

SimulationMSEl_bearing  =  0; 

SimulationMSE2_bearing  =  0; 

SimulationMSE3_bearing  =  0; 

for  k  =  1  :  L,  %  Monte  Carlo 

Simulation . 

rdNse  =  sqrt (nsePwr/2 )  *  randn (M, 1); 

rdNsel  =  sqrt (nsePwr/2 ) ; 

rd  =  rdo  +  rdNse (2 : end) -rdNse ( 1 ) ;  %  Noisy  source 

TDOAs  (range  differences) . 

ul  =  TDOALoc_LINEAR (S, rd, Q) ; 

[u2,  senor_pos]  = 

Linear_second_stage_polar_RDOA (atan2 (ul (2) , ul (1) ) , LM, LL/2, M 
LE_runs, [norm (uo, 2) , atan2 (uo (2) , uo (1) ) ] ' , rdNsel, [norm (ul, 2) 

, atan2  (ul  (2) ,ul (1) ) ]  ' ) ; 

u2  =  [u2 ( 1 ) *cos (u2  (2) ) ;u2  (1) *sin  (u2  (2) ) ] ; 

[u3,  senor_posl]  = 

Linear_second_stage_polar_RDOA (atan2 (ul (2) , ul (1) ) , LM, LL/2, M 
LE_runs+3 , [norm (uo, 2 ) , atan2 (uo (2 ) , uo ( 1 ) ) ] ' , rdNsel , [norm (ul , 
2) , atan2 (ul (2) , ul (1) ) ] ' ) ; 

u3  =  [ u  3 (1) *cos ( u  3 (2) ) ;u3 (1) *  s  i  n  ( u  3 (2) ) ] ; 

SimulationMSEl  =  SimulationMSEl  +  norm (ul-uo, 2 ) A2 ; 
SimulationMSE2  =  Simulat ionMSE2  +  norm (u2-uo, 2 ) A2 ; 
SimulationMSE3  =  SimulationMSE3  +  norm (u3-uo, 2 ) A2 ; 
SimulationMSEl_bearing  =  SimulationMSEl_bearing  + 
(atan2  (ul (2) , ul (1) ) -theta_tx) A2 ; 

SimulationMSE2_bearing  =  SimulationMSE2_bearing  + 
(atan2  (u2 (2) ,u2 (1) ) -theta  tx) A2 ; 
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SimulationMSE3_bearing  =  SimulationMSE3_bearing  + 
(atan2 (u3 (2) ,u3 (1) ) -theta_tx) A2 ; 
end; 

msel (nseldx)  =  sqrt (SimulationMSEl/L) ; 
mse2 (nseldx)  =  sqrt (SimulationMSE2/L) ; 
mse3 (nseldx)  =  sqrt (SimulationMSE3/L) ; 

msel_bearing (nseldx)  =  sqrt (SimulationMSEl_bearing/L) ; 
mse2_bearing (nseldx)  =  sqrt (SimulationMSE2_bearing/L) ; 
mse3_bearing (nseldx)  =  sqrt (SimulationMSE3_bearing/L) ; 

end; 

fprintf  (  ' \n ' ) ; 

figure (1) ;  plot (NsePwrVecdB, (msel ) , ' xk ' , ' Markers ize ' , 8 ) ; 
hold  on; 

plot (NsePwrVecdB, (mse2 ),' *r ',' MarkerSize 8 ) ;  hold  on; 
plot (NsePwrVecdB,  (crib) ,  '  k '  )  ; 
hold  on; 

plot (NsePwrVecdB, (crlb_LM) , ' --g ' ) ; 
grid  on;  hold  off; 

xlabel('RDOA  Noise  STD  dB (c\sigma) ' ) ; 
ylabel ( 'RMSE  Position  (m) ' ) ; 

legend ('1st  Phase  weighted  least-squares  Estimate 2nd 
Phase  ML  Estimate 1st  Phase  CRLB  RMSE  Position',  '2nd 
Phase  CRLB  RMSE  Position'); 

figure ( 3 ) ; 

plot (NsePwrVecdB, (msel_bearing) , ' xk ' , ' MarkerSize ' , 8) ;  hold 
on; 

plot (NsePwrVecdB, (mse2_bearing) , ' *r ' , ' MarkerSize ' , 8) ;  hold 
on; 

plot (NsePwrVecdB, (crlb_LM_bearing) , '-g'); 
grid  on;  hold  off; 

xlabel('RDOA  Noise  STD  dB (c\sigma) ' ) ; 
ylabel ( 'RMSE  Bearing  (rad) ' ) ; 

legend ('1st  Phase  weighted  least-squares  Estimate ',' 2nd 
Phase  ML  Estimate ',' 2nd  Phase  CRLB  RMSE  Bearing'); 
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Script  Name:  Signal  estimator  versus  wavelength  for  uniform  position  errors 

clc 

clear  all 
close  all 
num_run  =  1000; 
rand ( ' seed ' ,  0 ) 
randn ( ' seed '  ,  0 ) 
c  =  3e8; 

f  =  50:100:2600; 
f  =  f  *  le6; 
lambda  =  c  .  /  ( f )  ; 
k  =  2*pi . /lambda ; 

Number_of_real_nodes_Center  =  10; 

Number_of_real_nodes_Outer  =  Number_of_real__nodes_Center 
ML  =  Number_of_real_nodes_Center*3 ; 

%Signal  1 
Vs_l_real  =  1; 

Phase_l_real  =  degtorad(O); 

DOA_l_real  =  0; 

deploy_radius  =  200; 

X_center  =  0; 

Y_center  =  0; 

X_outer  =  deploy_radius ; 

Y_outer  =  0; 

X  =  0; 

Y  =  0; 

SNR  =  5; 

%%creating  random  receive  matrix 
a  =  -1/10; 
b  =  -a; 
aa  =  - . 9e-3 ; 
bb  =  -aa; 

Signal_estimates  =  zeros  (6,1); 

Signal_estimates_f rom_mean  =  zeros  (2,1); 
deltaX  =  b; 
deltaY  =  b; 

DOAn  =  bb ; 

Tans_weighted  least-squaresE  =  zeros  (2,1); 

I_Q_data  =  zeros (3,1); 

DOA  ESTIMATES  =  zeros (num  run, 1) ; 
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RMSE  =  zeros ( 4 ,  length (Number_of_real_nodes_Center )) ; 
num_pts  =  Number_of_real_nodes_Center ; 
for  MCR  =  1:  length (k) 


MCR 

%  this  is  the  -15  dB  run 
Tans_weighted  least-squaresE  =  0; 
for  xx  =  1 : num  run 

XnC  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

YnC  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

XnOL  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

YnOL  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

XnOR  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

YnOR  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

DOA_l_est  =  Locator_Function ( [ 2000 ; degtorad  ( 7 0 ) ] , - 
10)  -  degtorad ( 70 ) ; 

Rx_sig_center  =  Vs_l_real . *exp ( ( li ) . * (Phase_l_real 
+  k (MCR) . * ( (X_center  +  XnC) . *sin (DOA_l_est )  +  (Y_center  + 
YnC) . *cos (DOA_l_est) ) ) ) ; 

Rx_sig_center_w_AWGN  =  awgn (Rx_sig_center ,  SNR) ; 
%Summend  Center  Signal 
Center_IQ_Est  =  sum (Rx_sig_center ) ; 

Rx_sig_outer_lef t  = 

Vs_l_real . *exp ( (li) . * (Phase_l_real  +  k(MCR).*((- 
deploy_radius+  XnOL) . *sin (DOA_l_est )  +  (Y_outer  + 

YnOL) . *cos (DOA_l_est) ) ) ) ; 

Rx_sig_outer_lef t_w_AWGN  =  awgn (Rx_sig_outer_lef t , 

SNR)  ; 

%Estimate  Mag  and  Phase  for  Left  Cluster 
Left_IQ_Est  =  sum (Rx_sig_outer_lef t ) ; 

Rx_sig_outer_right  = 

Vs_l_real . *exp ( (li) . * (Phase_l_real  + 

k (MCR) . * ( (deploy_radius+  XnOR) . *sin (DOA_l_est )  +  (Y_outer  + 
YnOR)  . *  cos (DOA_l~est) ))); 

Rx_sig_outer_right_w_AWGN  = 
awgn (Rx_sig_outer_right ,  SNR) ; 

%Estimate  Mag  and  Phase  for  Right  Cluster 
Right_IQ_Est  =  sum (Rx_sig_outer_right ) ; 
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Sig_IQ_est  = 

abs (sum( [Center_IQ_Est ; Lef t_IQ_Est ; Right_IQ_Est ] ) ) /ML; 

Tans_weighted  least-squares  =  [Tans_weighted  least- 
squaresE  Sig_IQ_est] ; 

end 

Tans_weighted  least-squares  =  Tans_weighted  least- 
squaresE (2 : end) ; 

RMSE(1,MCR)  =  mean (Tans_weighted  least-squaresE) ; 

Tans_weighted  least-squares  =  0; 
for  xx  =  1 : num  run 

XnC  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

YnC  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

XnOL  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

YnOL  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

XnOR  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

YnOR  =  (a  +  (b-a) . *rand (num_pts , 1 ) ) ; 

DOA_l_est  =  Locator_Function ( [2000 ; degtorad (70 )], - 
5)  -  degtorad ( 7 0 ) ; 

Rx_sig_center  =  Vs_l_real . *exp ( ( li ) . * (Phase_l_real 
+  k (MCR) . * ( (X_center  +  XnC) . *sin (DOA_l_est )  +  (Y_center  + 
YnC) . *cos (DOA_l_est) ) ) ) ; 

Rx_sig_center_w_AWGN  =  awgn (Rx_sig_center ,  SNR) ; 
%Summend  Center  Signal 
Center_IQ_Est  =  sum (Rx_sig_center ) ; 

Rx_sig_outer_lef t  = 

Vs_l_real . *exp ( (li) . * (Phase_l_real  +  k (MCR) . * ( (- 
deploy_radius+  XnOL) . *sin (DOA_l_est )  +  (Y_outer  + 

YnOL) . *cos (DOA_l_est) ) ) ) ; 

Rx_sig_outer_lef t_w_AWGN  =  awgn (Rx_sig_outer_lef t , 

SNR)  ; 

%Estimate  Mag  and  Phase  for  Left  Cluster 
Left_IQ_Est  =  sum (Rx_sig_outer_lef t ) ; 

Rx_sig_outer_right  = 

Vs_l_real . *exp ( (li) . * (Phase_l_real  + 

k (MCR) . * ( (deploy_radius+  XnOR) . *sin (DOA_l_est )  +  (Y_outer  + 
YnOR)  . *  cos (DOA_l~est) ))); 

Rx_sig_outer_right_w_AWGN  = 
awgn (Rx_sig_outer_right ,  SNR) ; 
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%Estimate  Mag  and  Phase  for  Right  Cluster 
Right_IQ_Est  =  sum (Rx_sig_outer_right ) ; 

Sig_IQ_est  = 

abs (sum( [Center_IQ_Est ; Lef t_IQ_Est ; Right_IQ_Est ] ) ) /ML; 

Tans_weighted  least-squaresE  =  [Tans_weighted 
least-squaresE  Sig_IQ_est] ; 

end 

Tans_weighted  least-squaresE  =  Tans_weighted  least- 
squaresE  (2  :  end)  ; 

RMSE(2,MCR)  =  mean (Tans_weighted  least-squaresE); 

end 

figure ( 1 ) 

plot (lambda, RMSE (1, : ) , lambda, RMSE (2, : ) ) 

xlabel ( ' Lambda ( \lambda) ' ) ; 

ylabel ( 'Average  Normalized  Gain'); 

legend ('RDOA  noise  STD  (-10  dB)','RDOA  noise  STD  (-5 
dB)  '  )  ; 

figure (2 ) 

plot (f . /le6,RMSE (1, : ) , f . /le6,RMSE (2, : ) ) 

xlabel ( ' Lambda ( \ lambda) ' ) ; 

ylabel (' Average  Normalized  Gain'); 

legend ( ' RDOA  noise  STD  (-10  dB)','RDOA  noise  STD  (-5 
dB)  '  )  ; 
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